Foreword

This R Markdown document is a final output for Tree Phenology Analysis using R (Winter Semester 2022). I am completing this file while sitting on my hard chair during a cold night in March. This document will be a “documentation” of how I learned R for the first time truly. With the help of the course modules, notes from the lectures, and of Mr. Google.

I hope this work will be of use somehow to the course’s creators– a humble perspective from an R beginner who strove to learn things from bottom upwards. For sure, I will be forever grateful to the existence of this course because it introduced me to the tool I hope I will be using often in my career. Prost!

Additionally, since I am a self-proclaimed nerd who likes thoroughly dissecting concepts, some chapters will start with condensed descriptions of the themes discussed in class (and I use mathematical and logic operation symbols often to show relationships), but these are not complete regurgitation of what was discussed; some are seemingly trivial terms which are new to me, followed by the important bit– the exercises and my take on them.

I hope this work of an R newbie is viewed with kindness and consideration. Admittedly, most of the time I just followed step by step the templates given as examples, and was not able to play around much with R features yet. As I came across problems, I attempted to troubleshoot and improvise what works and what does not, in order to progress in the exercises. My main aim in taking this course is to learn how R works in the context of scientific research, and in the end I think I achieved it somehow and more.

03 Tree dormancy

Dormancy prevents injury during winter. Triggered by environmental signals [photo-periods f(latitudes) + temperature f(geographical location, photoperiods, climate)], works through interaction of physiological processes.

Stages of dormancy:

  1. Dormancy establishment- environment signals (declining temp and photoperiod)

  2. Endo dormancy- tree’s endogenous properties dictate this (chilling temp)

  3. Eco dormancy-acclimatization but not deeply dormant, no morphological /development changes. temperature- main environmental driver

    (warm temperatures)

  4. growth resumption

Physiological processes regulating dormancy:

Dormancy determination

Dormancy-> prevents injuries. Cold is required to regain the capabilities to grow

Endo-dormancy, eco-dormancy (growth is prevented in unfavorable condition, heat exposure)

Empirical-> growth chambers, shoots with buds collected in winter. 7 to 12 days. Buds growth or not. Longer exposure to chilling? Chill period, forcing period

Statistical-> Long phenological datasets ) temp records. Based of long series of dates, previous temperature records. Calculates max, mean temperatures, flowering dates. Chilling and heat accumulation. Flowering dates with mean and max temperatures. Number of independent variables> dependent variables (for PLS)

Chilling period-> when model coefficients are + significant values (temperature delays bloom dates)

Forcing period-> when model coefficients are -significant values (less temperature delays bloom dates)

experimental (empirical)

statistical- (dormancy determination) research purposes, collected in research centers. new commercial cultivars. Endodormancy break, chilling hours

Phenological datasets

Phenology in practical agriculture-> irrigation scheduling (when water is mostly taken up)

Tree-> budding, flowering, fruiting stages

Old vs BBCH scales

  1. letters from A to J- every 2 days during spring. earliest and latest state

  2. limitations; need for common phenology across the world

    -need for standardized scale

    -use of numbers is more practical than letters

Scales for plant growth stages

Cereals (BBCH, monocot, dicot)

Exercise:

  1. Put yourself in the place of a breeder who wants to calculate the temperature requirements of a newly released cultivar. Which method will you use to calculate the chilling and forcing periods? Please justify your answer.

    To calculate temperature requirements of a newly-released cultivar, I would do empirical first (trial and error, actual field experiments); but this could take time and could be expensive because the researcher would have to wait seasons by seasons and replicate to have enough repetitions, simulate the conditions through controlled environment (growth chamber). Statistical can be done with available historical data (long phenological data), and it has to be “applicable” to the new cultivar, meaning, more similarities than differences. However, the two methods could be done in parallel since they require different sets of data anyway. Results can be compared to verify if the models generated through statistical method is acceptable enough.

  2. Which are the advantages (2) of the BBCH scale compared with earlier scales?

    a. Standardized across the globe

    b. More specific (has primary and secondary classifications in contrast to old version that used single letters)

  3. Classify the following phenological stages of sweet cherry according to the BBCH scale:

    budding- 54

    flowering- 65

    fruiting- 87

04 Climate change and impact projection

Drivers of climate change

these factors are quantifiable: attribution studies

models have been crafted, compared to observed, and the recent increase(since industrialization) in global temperature can be explained by GHGE from human activities (burning of fossil fuels)

trends of CO2 emissions is still going up (new climatic territories)

degree of CO2 concentrations never seen by mankind, at high rates

Recent warming

Starting around 1980s, there is dramatic rise in global temperatures. Globally warmest years are very recent.

-> Siberia (2020) >8C warmer (possibility of +feedback=acceleration of warming effect)

Maunder minimum-> minimal sunspots, low solar radiation

Industrialization- GHGEs

Global surface temperatures started with warming, then a long period of cooling.

+1.5C is optimistic expectation (global temperature increase), but still requires adaptation-> +5C would cause major global changes unseen before by mankind

Rainfall is harder to predict than temperatures (no clear trends)

Moisture deficit in EU

Places with airports= good data (weather stations)

Future Scenarios

No way of validating which scenario is correct

Future temperatures

RCP 8.5: rise of 2 to 11C in global temperatures

Tipping point- point of sudden transition to new state (poorly considered in impact projections)

-> meltwater from Greenland dilutes salt concentration, less driver for global ocean currents (cold, salty water is heaviest); further slowing= transition of global climate to new state

-Greenland, arctic iceloss-> weakened albedo

-Permafrost thawing-> CO2 and methane release (more warming)

full of uncertainties

emissions can be translated to temperature increase

Paris Agreement: 1.5C increase only

-> 2060, CO2 neutral atleast

Impact Projection

system response

Common limitations: climate data are scarce or poor quality; non-climatic factors also change through time, and difficult to project; CO2 can only be included in process-based models

Projecting climate-change impacts in complex agricultural systems:

need to understand what is happening: weather, climate, and crops + all other relevant factors affecting system performance; consider and communicate uncertainties

Exercise:

  1. Main drivers of climate change (decade to century scale), explain mechanism through which the currently most important driver affects our climate.

    at century scale: Sun (radiation), Aerosols (liquid,solid, mixed suspensions), Clouds (cooling and warming effect), Ozone (destroyed by CFCs, Montreal protocol of 1987), Surface albedo (reflects radiation) , GHGs, long term drivers (solar activity trends, ocean currents, plants vs animals, volcanic and meteorite activity, milankovic cycles)

    mechanism through which most important driver affects climate: rising concentration of GHGs-> this is due to activities of mankind. CO2, CH4, and SO2 are released into the atmosphere. These stay in the atmosphere for very long periods and absorb long-wave radiation, which causes heating up of the Earth; thus causing “greenhouse effect”.

  2. Explain briefly what is special about temperature dynamics of the recent decades, and why we have good reasons to be concerned.

    Industrial revolution led to sharp rise in global temperatures (land, ocean surface temperatures) never seen before by mankind. A world that is quickly becoming >1.5C hotter than now would require tremendous adaption effort from future generations. Since we have not seen it before, the fear of uncertainty and looming feeling of doom is a great motivation to take GHGE reduction seriously; or else we might say ‘hi’ to dinosaurs sooner than we think.

  3. What does the abbreviation ‘RCP’ stand for, how are RCPs defined, and what is their role in projecting future climates?

    RCP: Representative Concentration Pathways. Each RCP is defined by additional radiative forcing; the higher the number, the higher the emissions and global heating. Currently, the commonly used ones are: RCP 2.6, 4.5, 6.0, 8.5, order is from the most conservative to the most emission-intensive estimate.

    Impact projection starts at identifying RCPs -> inputs for climate change impact models. Process is as follows:

    RCPs> GCM > downscaling (RCM) > statistical downscaling>Temperature, Precipitation, CO2 projections > Impact projection

  4. Briefly describe the 4 climate impact projection methods described in the fourth video.

    a. Statistical models- climate parameters x impact measure; used to explain past trends and project future impacts

    b. Species distribution modelling- climate parameters x presence/absence of species.

    c. Process-based models- use of mathematical equations to quantify system components and processes involved, and to project climate impact to the system. Complex since all relevant factors must be included (non-climatic); resource-intensive and requires processing of huge data

    d. Climate analogue models- involves the use of a similar and existing location where the climate projection can be compared to. Target x analog locations

05 Winter chill projections

Chilling hours- hours of low temperature required to break fruit dormancy during winter. Various models have different definitions (Kaufmann, 2019):

Apples- 500 to 1000 CH

1 CP= 10 CH

Exercise:

  1. Sketch out three data access and processing challenges that had to be overcome in order to produce chill projections with state-of-the-art methodology.

    -limitation of equipment/available infrastructure in processing huge chunk of data

    -incomplete or lack of high resolution climate datasets for some places; would need to make estimates (ground for errors) with respect to nearby places with available data

    -form in which data is available in (especially for really large data); data in raster form adds additional steps before they can be rendered usable

  2. Outline, in your understanding, the basic steps that are necessary to make such projections.

    understand context and systems in location chosen (research + literature review)-> form your research question -> choose the appropriate chill model and climate projection ensembles to use -> obtain data needed for the chosen model-> interpolate/make estimates if there gaps -> process data using the model -> compare projections with existing projections/historic observations -> conclusion + expression of uncertainties and errors

06 Manual chill analysis

CH computation requirement: hourly data;

read.table or read.csv functions;

.xls or .xlsx winters_hours_gaps dataset kable format for improved aesthetic;

<- assigning new values to R objects comparison commands (<, <=, ==, => and >);

c() command creates strings of numbers “vectors”;

& command combines comparisons;

FALSE and TRUE ~ 0 and 1 which command to identify items in dataset of interest

Creating functions:

OurFunctionName <- function(argument1, argument2, ...) {ourCode}

Exercises:

  1. Write a basic function that calculates warm hours (>25°C)

    Year<-c(2000:2010)
    Temps<-c(20:30)
    randomdataset<-data.frame(Year, Temps)
    
    ##function creation, assigning of threshold to variable
    WH<- function(randomdataset)
    {
      threshold_warm<- 25
      randomdataset[,"Warm_year"]<- randomdataset$Temps>threshold_warm
      return(randomdataset)
    
    }
    
    sampleWH<-WH(randomdataset)
    
    write.csv(sampleWH, "data/sampleWH.csv", row.names = FALSE)

    Year

    Temps

    Warm_year

    2000

    20

    FALSE

    2001

    21

    FALSE

    2002

    22

    FALSE

    2003

    23

    FALSE

    2004

    24

    FALSE

    2005

    25

    FALSE

  2. Apply this function to the Winters_hours_gaps dataset

    library(knitr)
    library(pander)
    library(kableExtra)
    library(chillR)
    
    ##function creation applied to Winter_hours_gaps; I made a new column that tells whether the temperatures are greater than the set threshold value: 25C
    hourtemps<- Winters_hours_gaps[,c("Year", "Month", "Day", "Hour", "Temp")]
    
    WH<- function(hourtemps)
    {
      threshold_warm<- 25
      hourtemps[,"Warm_Hour"]<- hourtemps$Temp>threshold_warm
      return(hourtemps)
    }
    
    write.csv(WH(hourtemps),"data/WH.csv", row.names = FALSE)

    Year

    Month

    Day

    Hour

    Temp

    Warm_Hour

    2008

    3

    3

    10

    15.127

    FALSE

    2008

    3

    3

    11

    17.153

    FALSE

    2008

    3

    3

    12

    18.699

    FALSE

    2008

    3

    3

    13

    18.699

    FALSE

    2008

    3

    3

    14

    18.842

    FALSE

    2008

    3

    3

    15

    19.508

    FALSE

  3. Extend this function, so that it can take start and end dates as inputs and sums up warm hours between these dates

       library(chillR)
    
     ##function creation applied to Winter_hours_gaps; I made a new column that tells whether the temperatures are greater than the set threshold value: 25C
    hourtemps<- Winters_hours_gaps[,c("Year", "Month", "Day", "Hour", "Temp")]
    
    WH<- function(hourtemps)
    {
      threshold_warm<- 25
      hourtemps[,"Warm_Hour"]<- hourtemps$Temp>threshold_warm
      return(hourtemps)
    }
    
    WH(hourtemps)[13:20,]
    
    
    ##sum of warm hours of start date and end year in YEARMODA
    
    sum_WH<-function(hourtemps, Start_YEARMODA, End_YEARMODA)
    {
       Start_Year<-trunc(Start_YEARMODA/10000)
      Start_Month<- trunc((Start_YEARMODA-Start_Year*10000)/100)
      Start_Day<- Start_YEARMODA-Start_Year*10000-Start_Month*100
      Start_Hour<-12
      End_Year<-trunc(End_YEARMODA/10000)
      End_Month<- trunc((End_YEARMODA-End_Year*10000)/100)
      End_Day<- End_YEARMODA-End_Year*10000-End_Month*100
      End_Hour<-12
    
      Start_Date<-which(hourtemps$Year==Start_Year & hourtemps$Month==Start_Month & hourtemps$Day==Start_Day & hourtemps$Hour==Start_Hour)
      End_Date<-which(hourtemps$Year==End_Year & hourtemps$Month==End_Month & hourtemps$Day==End_Day & hourtemps$Hour==End_Hour)
    
      Warm_hours<- WH(hourtemps)
      return(sum(Warm_hours$Warm_Hour[Start_Date:End_Date]))
    }
    ##sample content for checking. # of warm hours in set time interval by manual counting= using sum_WH function), so I confirmed there's no error anywhere
    sum_WH(hourtemps,20080517,20080518)
    WH(hourtemps)[1803:1827,]

    Year

    Month

    Day

    Hour

    Temp

    Warm_Hour

    1803

    2008

    5

    17

    12

    31.663

    TRUE

    1804

    2008

    5

    17

    13

    32.639

    TRUE

    1805

    2008

    5

    17

    14

    33.183

    TRUE

    1806

    2008

    5

    17

    15

    33.417

    TRUE

    1807

    2008

    5

    17

    16

    31.791

    TRUE

    1808

    2008

    5

    17

    17

    29.790

    TRUE

    1809

    2008

    5

    17

    18

    29.790

    TRUE

    1810

    2008

    5

    17

    19

    27.087

    TRUE

    1811

    2008

    5

    17

    20

    28.196

    TRUE

    1812

    2008

    5

    17

    21

    28.493

    TRUE

    1813

    2008

    5

    17

    22

    27.875

    TRUE

    1814

    2008

    5

    17

    23

    21.843

    FALSE

    1815

    2008

    5

    18

    0

    19.175

    FALSE

    1816

    2008

    5

    18

    1

    18.176

    FALSE

    1817

    2008

    5

    18

    2

    18.533

    FALSE

    1818

    2008

    5

    18

    3

    16.201

    FALSE

    1819

    2008

    5

    18

    4

    15.008

    FALSE

    1820

    2008

    5

    18

    5

    14.721

    FALSE

    1821

    2008

    5

    18

    6

    15.199

    FALSE

    1822

    2008

    5

    18

    7

    17.938

    FALSE

    1823

    2008

    5

    18

    8

    23.088

    FALSE

    1824

    2008

    5

    18

    9

    25.162

    TRUE

    1825

    2008

    5

    18

    10

    27.161

    TRUE

    1826

    2008

    5

    18

    11

    29.615

    TRUE

    1827

    2008

    5

    18

    12

    31.459

    TRUE

    ##[1] 15

07 Chill models

CH model- not credible

Utah Model- more credible (different weights for different temperatures

Dynamic Model- complex, needs heavy maths.

Chilling_Hours() function determines if each hour are within 0 to 7.2C ;

step_model() function- allows user to define model (set thresholds and weights);

summ==TRUE cumulative sum;

data.frame- input, contains weights for temperature ranges; custom() function implements chill model based on set;

data.frame input make_JDay() function adds Julian dates;

chilling() function implements CH, Utah, and Dynamic models, and computes GDH (Growing Degree Days);

tempResponse function lets user pick what model they want in chilling() function

Exercise:

  1. Run the chilling() function on the Winters_hours_gap dataset

    library(chillR)
    chill_sample<-chilling(make_JDay(Winters_hours_gaps),Start_JDay = 90, End_JDay = 100)
    
    write.csv(chill_sample, "data/chill_sample.csv", row.names = FALSE)

    Season

    End_year

    Season_days

    Data_days

    Perc_complete

    Chilling_Hours

    Utah_Model

    Chill_portions

    GDH

    2007/2008

    2008

    11

    11

    100

    40

    15.5

    2.009147

    2406.52

  2. Create your own temperature-weighting chill model using the step_model() function

    I looked at Milech et. al (2018) to see other existing models I have not seen yet so far in this module, and saw a nice table from their study summarizing models and the ranges. I picked Low Chill Model to be my guinea pig for this exercise portion; I did not use its equations, only the temperature ranges and weights.

    ##implementing the ranges from Low Chill Model
    
    df<-data.frame(
      lower= c(-1000, 1.8, 7.9, 13.9, 16.9, 19.4, 20.4),
      upper=c(1.8, 7.9,13.9,16.9, 19.4, 20.4, 1000),
      weight=c(0, 0.5, 1, 0.5, 0, -0.5, -1)
    )
    
    custom<- function(x) step_model(x, df)
    
    ##sample values for first 100 rows from the Winters_hours_gaps dataset)
    custom(Winters_hours_gaps$Temp)[1:100]
    
    write.csv(custom(Winters_hours_gaps$Temp), "data/custom.csv", row.names = FALSE)
  3. Run this model on the Winters_hours_gaps dataset using the tempResponse() function.

    df<-data.frame(
      lower= c(-1000, 1.8, 7.9, 13.9, 16.9, 19.4, 20.4),
      upper=c(1.8, 7.9,13.9,16.9, 19.4, 20.4, 1000),
      weight=c(0,0.5, 1, 0.5, 0, -0.5, -1)
    )
    
    custom<- function(x) step_model(x, df)
    
    output<-tempResponse(make_JDay(Winters_hours_gaps), Start_JDay = 90, End_JDay = 100, models=list(custom=custom))
    
    
    write.csv(output, "data/output.csv", row.names = FALSE)

    Season

    End_year

    Season_days

    Data_days

    Perc_complete

    custom

    2007/2008

    2008

    11

    11

    100

    96.5

08 Making hourly temperatures

daylength function- generates daily temperature curves from the latitude input;

ggplot2 package- for plotting;

JDay= Julian Date (1 to 365 or 366 in leap years);

melt command of the reshape2 package- separates/melts dataframe, stacked; stack_hourly_temps() function integrates daily dynamics using minimum and maximum temperature data + latitude;

KA_weather dataset contains temperature data from Bonn Uni; Empirical_daily_temperature_curve() function determines typical temperatures at locations; make_all_day_table function fills gaps in daily and hourly temperature records

Exercise:

  1. Choose a location of interest, find out its latitude and produce plots of daily sunrise, sunset and daylength.

    I have always been fascinated with Nordic mythology– Vikings. Sometimes I ask myself: why don’t we just consult the three Norns who spin the webs of fate: Urd (Past), Verdandi (Present), and Skuld (future) instead of modelling climate scenarios? Would be a cooler way to predicting the heating of Earth! Below is a fun image of the three Norns I got from a Norse Mythology website.

    My Nordic tale-loving soul convinced me to look at interesting locations in Norway, Sweden, and Denmark. I narrowed down my search by looking at the ranking of countries with the most number of weather stations, and voila! got my answer: Sweden. The plot below shows the number of weather stations across the globe from Jaffre (2019).

    An interesting project implemented in Northern Sweden by Brännland Cider (a Swedish company making ice cider from apples) caught my attention. European innovation partnership for productivity and sustainability within agriculture (EIP-agri) of the Swedish Board of Agriculture granted them around 935,000 euro for this. The aim is to pilot climate-resilient agriculture in the Northern EU by trying to grow apples commercially in the harsh northern climate. The three counties involved are: Norrbotten, Västerbotten and Västernorrland; and I picked Västerbotten (64.34° N, 18.31° E) because it had many airports close to it– nearest was 12 km away. It would be very interesting to know what the field trials would look years in the future temperature-wise. Would they be able to make ice cider decades from now in that area?

    I kept in mind that the location will be used repeatedly for the next topics so I chose quite strategically. Though my expectation is that they will be able to continue to grow apples because that region will not get warmer as per previous projections by Luedeling et al. (2011):

    Now, for the exercise,

    require(chillR)
    require(ggplot2)
    require(reshape2)
    Days<-daylength(latitude=64.34, JDay=1:365)
    Days_df<-data.frame(JDay=1:365, Sunrise=Days$Sunrise, Sunset=Days$Sunset, Daylength=Days$Daylength)
    Days_df<-melt(Days_df, id=c("JDay"))
    ggplot(Days_df, aes(JDay, value))+ geom_line(lwd=1.5) + facet_grid(cols=vars(variable))+ylab("Time of Day/Daylength (Hours)")+ theme_bw(base_size=20)

  2. Produce an hourly dataset, based on idealized daily curves, for the KA_weather dataset (included in chillR)

    require(chillR)
    
    KA_weather_stacked<-stack_hourly_temps(KA_weather, latitude=64.34)
    
    write.csv(KA_weather_stacked, "data/KA_weather_stacked.csv", row.names = FALSE)

    hourtemps.Year

    hourtemps.Month

    hourtemps.Day

    hourtemps.Tmax

    hourtemps.Tmin

    hourtemps.JDay

    hourtemps.Hour

    hourtemps.Temp

    QC

    1998

    1

    1

    8.2

    5.1

    1

    0

    5.1

    NA

    1998

    1

    1

    8.2

    5.1

    1

    1

    5.1

    NA

    1998

    1

    1

    8.2

    5.1

    1

    2

    5.1

    NA

    1998

    1

    1

    8.2

    5.1

    1

    3

    5.1

    NA

    1998

    1

    1

    8.2

    5.1

    1

    4

    5.1

    NA

    1998

    1

    1

    8.2

    5.1

    1

    5

    5.1

    NA

  3. Produce empirical temperature curve parameters for the Winters_hours_gaps dataset, and use them to predict hourly values from daily temperatures (this is very similar to the example above, but please make sure you understand what’s going on)

    library(chillR)
    library(ggplot2)
    
    ##determining empirical daily temperatures from hourly data using Winters_hours_gaps dataset
    coeffs<-Empirical_daily_temperature_curve(Winters_hours_gaps)
    Winters_daily<-make_all_day_table(Winters_hours_gaps, input_timestep="hour")
    Winters_hours<-Empirical_hourly_temperatures(Winters_daily,coeffs)
    
    ?Winters_hours_gaps
    
    ##input latitude of Winters, California (38.5 N) where Winters_hours_gaps dataset was recorded in
    require(reshape2)
    
    Winters_hours<-Winters_hours[,c("Year","Month","Day","Hour","Temp")]
    colnames(Winters_hours)[ncol(Winters_hours)]<-"Temp_empirical"
    Winters_ideal<-stack_hourly_temps(Winters_daily, latitude=38.5)$hourtemps
    Winters_ideal<-Winters_ideal[,c("Year","Month","Day","Hour","Temp")]
    colnames(Winters_ideal)[ncol(Winters_ideal)]<-"Temp_ideal"
    
    ##Triangular dataset
    Winters_triangle<-Winters_daily
    Winters_triangle[,"Hour"]<-0
    Winters_triangle$Hour[nrow(Winters_triangle)]<-23
    Winters_triangle[,"Temp"]<-0
    Winters_triangle<-make_all_day_table(Winters_triangle,timestep="hour")
    colnames(Winters_triangle)[ncol(Winters_triangle)]<-"Temp_triangular"
    
    for(i in 2:nrow(Winters_triangle))
    {if(is.na(Winters_triangle$Tmin[i])) Winters_triangle$Tmin[i]<-Winters_triangle$Tmin[i-1]
    if(is.na(Winters_triangle$Tmax[i])) Winters_triangle$Tmax[i]<-Winters_triangle$Tmax[i-1]
    }
    Winters_triangle$Temp_triangular<-NA
    
    Winters_triangle$Temp_triangular[which(Winters_triangle$Hour==6)]<-
      Winters_triangle$Tmin[which(Winters_triangle$Hour==6)] 
    Winters_triangle$Temp_triangular[which(Winters_triangle$Hour==18)]<-
      Winters_triangle$Tmax[which(Winters_triangle$Hour==18)] 
    Winters_triangle$Temp_triangular<-
      interpolate_gaps(Winters_triangle$Temp_triangular)$interp
    Winters_triangle<-Winters_triangle[,c("Year","Month","Day","Hour","Temp_triangular")]
    
    ##merging the dataframes
    Winters_temps<-merge(Winters_hours_gaps,Winters_hours, by=c("Year","Month","Day","Hour"))
    Winters_temps<-merge(Winters_temps,Winters_triangle, by=c("Year","Month","Day","Hour"))
    Winters_temps<-merge(Winters_temps,Winters_ideal, by=c("Year","Month","Day","Hour"))
    
    ##conversion of dates, reorganizing dataframes, plotting
    Winters_temps[,"DATE"]<-ISOdate(Winters_temps$Year,Winters_temps$Month, Winters_temps$Day, Winters_temps$Hour)
    
    
    Winters_temps_to_plot<-Winters_temps[,c("DATE","Temp","Temp_empirical","Temp_triangular","Temp_ideal")]
    Winters_temps_to_plot<-Winters_temps_to_plot[100:200,]
    Winters_temps_to_plot<-melt(Winters_temps_to_plot, id=c("DATE")) 
    colnames(Winters_temps_to_plot)<-c("DATE","Method","Temperature")
    
    
    ggplot(data=Winters_temps_to_plot, aes(DATE,Temperature, colour=Method)) +
      geom_line(lwd=1.3) + ylab("Temperature (°C)") + xlab("Date")

09 Getting temperature data

GSOD- Global Summary of Day database

ChillR only has access to GSOD and California database

handle_gsod() function for data retrieval from GSOD;

Overlap_years column shows # years available;

Perc_interval_covered column shows % that target (start-end date input) was covered

Exercise:

  1. Choose a location of interest and find the 25 closest weather stations using the handle_gsod function

    Will use Västerbotten, Sweden (64.34° N, 18.31° E) again.

    library(chillR)
    ##generation of station list based on target location's coordinates
    station_list_Västerbotten<-handle_gsod(action="list_stations",
                              location=c(18.31,64.34),
                              time_interval=c(1990,2020))
    
    write.csv(station_list_Västerbotten, "data/station_list_Västerbotten.csv")

    chillR_code

    STATION.NAME

    CTRY

    Lat

    Long

    BEGIN

    END

    distance

    Overlap_years

    Perc_interval_covered

    022630_99999

    FREDRIKA

    SW

    64.083

    18.367

    20010701

    20230307

    28.78

    19.50

    62.91

    022640_99999

    FREDRIKA

    SW

    64.083

    18.417

    19730101

    20021028

    29.12

    12.83

    41.37

    023260_99999

    FREDRIKA

    SW

    64.083

    18.417

    19770701

    20041221

    29.12

    14.97

    48.30

    022455_99999

    LYCKSELE

    SW

    64.548

    18.716

    19980408

    20230307

    30.33

    22.73

    73.33

    022610_99999

    LYCKSELE

    SW

    64.550

    18.717

    19730101

    20230307

    30.53

    31.00

    100.00

    020630_99999

    ASELE

    SW

    64.167

    17.367

    19550101

    19770630

    49.62

    0.00

    0.00

  2. Download weather data for the most promising station on the list

    From the 25 stations nearest the chosen location, the fifth (Lycksele, ChillR code 022610-99999) was the most promising; it covered 100% of the data for 1990 to 2020. Though it was quite far in terms of distance from the target area.

    library(chillR)
    
    ##generation of station list based on target location's coordinates
    station_list<-handle_gsod(action="list_stations",
                              location=c(18.31,64.34),
                              time_interval=c(1990,2020))
    
    ##having to enter just the number of station from the list created was a relief, thanks chillR coders!
    weather<-handle_gsod(action="download_weather",
                         location=station_list$chillR_code[5],
                         time_interval=c(1990,2020))
    
    ##check sample data
    weather[[1]]$weather[5090:6010,]
    
    #downloading data
    write.csv(weather,"data/Västerbotten_weather.csv",row.names=FALSE)

    LYCKSELE.database

    LYCKSELE.weather.DATE

    LYCKSELE.weather.STN…

    LYCKSELE.weather.WBAN

    LYCKSELE.weather.YEAR

    LYCKSELE.weather.MONTH

    LYCKSELE.weather.DAY

    LYCKSELE.weather.TEMP

    LYCKSELE.weather.Count1

    LYCKSELE.weather.DEWP

    LYCKSELE.weather.Count2

    LYCKSELE.weather.SLP

    LYCKSELE.weather.Count3

    LYCKSELE.weather.STP

    LYCKSELE.weather.Count4

    LYCKSELE.weather.VISIB

    LYCKSELE.weather.Count5

    LYCKSELE.weather.WDSP

    LYCKSELE.weather.Count6

    LYCKSELE.weather.MXSPD

    LYCKSELE.weather.GUST

    LYCKSELE.weather.MAX

    LYCKSELE.weather.MaxFlag

    LYCKSELE.weather.MIN

    LYCKSELE.weather.MinFlag

    LYCKSELE.weather.PRCP

    LYCKSELE.weather.PrcpFlag

    LYCKSELE.weather.SNDP

    LYCKSELE.weather.FRSHTT

    LYCKSELE.weather.Year

    LYCKSELE.weather.Month

    LYCKSELE.weather.Day

    GSOD

    1990-01-01 12:00:00

    22610

    99999

    1990

    1

    1

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    1990

    1

    1

    GSOD

    1990-01-02 12:00:00

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    1990

    1

    2

    GSOD

    1990-01-03 12:00:00

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    1990

    1

    3

    GSOD

    1990-01-04 12:00:00

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    1990

    1

    4

    GSOD

    1990-01-05 12:00:00

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    1990

    1

    5

    GSOD

    1990-01-06 12:00:00

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    NA

    1990

    1

    6

  3. Convert the weather data into chillR format

    library(chillR)
    
    ##generation of station list based on target location's coordinates
    station_list<-handle_gsod(action="list_stations",
                              location=c(18.31,64.34),time_interval = c(1990,2020))
    
    ##having to enter just the number of station from the list created was a relief, thanks chillR coders!
    weather<-handle_gsod(action="download_weather",
                         location=station_list$chillR_code[5],
                         time_interval=c(1990,2020))
    
    
    ##datacleaning, changing format to chillR
    cleaned_weather<-handle_gsod(weather)
    
    ##check sample cleaned data, I used a non-boring range
    cleaned_weather[[1]]$weather[5090:6010,]
    
    
    ##downloading data for future use
    write.csv(cleaned_weather$weather,"data/Västerbotten_chillR_weather.csv",row.names=FALSE)

From the looks of the data I got as shown above from Västerbotten taken from GSOD database, it’s not possible to work with it for the next exercises. I already tried doing first exercise in Chapter 10 and it failed to be read using read.csv function (contains no data, that also happened to Bonn weather example from GSOD I ran). So for the next exercises, I would have to choose another location, possibly in Germany since the DWD database seems working well with chillR. To be able to move forward, I will also use DWD database for the filling of gaps exercise.

Another reason is we already know that Nordic regions will remain cold enough anyway as per previous projections. Makes me sad though; but will not cry because I am a tough Viking by heart.

10 Getting temperature data

linear interpolation can be done to fill in short gaps (hours to weeks), but longer gaps need better methods to reflect temperature dynamics

interpolate_gaps() function carries out simple interpolation; fix_weather() function uses linear interpolation between Tmax and Tmin columns and can be used to identify gaps in the dataset

Exercise:

  1. Use chillR functions to find out how many gaps you have in this dataset (even if you have none, please still follow all further steps)

    I will repeat the previous steps done in the last exercise in this portion. New location of choice is the Swabian Orchard Paradise (48.4931° N, 9.3998° E) near Stuttgart. I had the chance to visit a Schwabian place before and ate some of their local dishes, courtesy of the tour by my mentor Martin Gummert. The orchard grows cherries, pears, and plums– and it is said that it has been centuries-old; would be interesting to check if they will survive another century.

    A quick view of the place I got from Schwäbisches Streuobstparadies Website (n.d.):

    Now onto the code:

    library(chillR)
    
    ##generation of list, handle DWD dataset
    station_list_dwd<-handle_dwd(action="list_stations",
                              location=c(9.40,48.49),
                              time_interval=c(19900101,20201231))
    
    ##check weather station list
    require(kableExtra)
    
    kable(station_list_dwd) %>%
      kable_styling("striped", position = "left", font_size = 8)
    
    ##the list showed M<fc>nsingen-Apfelstetten (station 654) as the nearest and most number of years available (1946 to 2021), so I am choosing it (#2 in the station list)
    
    
    ##downloading the data from chosen weather station and checking data
    weather_dwd<-handle_dwd(action="download_weather",
                         location=station_list_dwd$Station_ID[2],
                         time_interval=c(19900101,20201231))
    
    head(weather_dwd$`M<fc>nsingen-Apfelstetten`$weather)
    
    ##datacleaning and checking data
    cleaned_weather_dwd<-handle_dwd(weather_dwd)
    
    cleaned_weather_dwd$`M<fc>nsingen-Apfelstetten`
    
    
    ##saving data in computer
    write.csv(station_list_dwd,"data/station_list_dwd.csv",row.names=FALSE)
    write.csv(weather_dwd,"data/MApfelstetten_weather1.csv",row.names=FALSE)
    write.csv(cleaned_weather_dwd$`M<fc>nsingen-Apfelstetten`,"data/MApfelstetten_chillR_weather.csv",row.names=FALSE)
    
    ##assigning previous data to an object 
    MApfelstetten<-read.csv("data/MApfelstetten_chillR_weather.csv")
    
    #checking for gaps using fix_weather function
    MApfelstetten_QC<-fix_weather(MApfelstetten)$QC
    
    write.csv(MApfelstetten_QC, "data/MApfelstetten_QC.csv", row.names = FALSE)
    
    
    ##saving and downloading the needed columns for next exercises
    MApfelstetten_weather<-MApfelstetten[,c("Year","Month", "Day", "Tmax", "Tmin")]
    kable(MApfelstetten_weather,) %>%
      kable_styling("striped", position = "left", font_size = 10)
    
    write.csv(MApfelstetten_weather, "data/MApfelstetten_weather.csv", row.names = FALSE)

    Quality control summary produced by fix_weather()

    Season

    End_year

    Season_days

    Data_days

    Missing_Tmin

    Missing_Tmax

    Incomplete_days

    Perc_complete

    1989/1990

    1990

    365

    365

    0

    0

    0

    100

    1990/1991

    1991

    365

    365

    0

    0

    0

    100

    1991/1992

    1992

    366

    366

    0

    0

    0

    100

    1992/1993

    1993

    365

    365

    0

    0

    0

    100

    1993/1994

    1994

    365

    365

    0

    0

    0

    100

    1994/1995

    1995

    365

    365

    0

    0

    0

    100

    Year

    Month

    Day

    Tmax

    Tmin

    1990

    1

    1

    0.6

    -10.0

    1990

    1

    2

    4.3

    -8.8

    1990

    1

    3

    -2.3

    -4.7

    1990

    1

    4

    -3.3

    -5.5

    1990

    1

    5

    0.4

    -8.5

    1990

    1

    6

    0.5

    -2.5

    Surprise, surprise! according to the results from running the fix_weather() function, there is no gap for the Münsingen-Apfelstetten weather station (DWD database) for our years of interest.

    So now, I will create artificial gaps by deleting some data from the original dataset and work with filling it for the next exercises. The dataset I am using is here.

    library(chillR)
    
    ##assigning previous data to an object 
    MApfelstetten_gaps<-read.csv("data/MApfelstetten_chillR_weather_gaps.csv")
    
    #checking for gaps using fix_weather function
    MApfelstetten_QC_gaps<-fix_weather(MApfelstetten_gaps)$QC
    
    write.csv(MApfelstetten_QC_gaps, "data/MApfelstetten_QC_gaps.csv", row.names = FALSE)

    Quality control summary produced by fix_weather()

    Season

    End_year

    Season_days

    Data_days

    Missing_Tmin

    Missing_Tmax

    Incomplete_days

    Perc_complete

    1989/1990

    1990

    365

    365

    38

    38

    38

    89.6

    1990/1991

    1991

    365

    365

    42

    42

    42

    88.5

    1991/1992

    1992

    366

    366

    0

    0

    0

    100.0

    1992/1993

    1993

    365

    365

    0

    0

    0

    100.0

    1993/1994

    1994

    365

    365

    0

    0

    0

    100.0

    1994/1995

    1995

    365

    365

    0

    0

    0

    100.0

    This confirmed the existence of gaps in 1990, 1991, 2004, 2005, 2012, 2018, 2019, and 2020.

  2. Create a list of the 25 closest weather stations using the handle_gsod function

    I will use the location of the Swabian Orchard Paradise (48.4931° N, 9.3998° E) to check weather stations around via the handle_dwd function.

    library(chillR)
    library(kableExtra)
    
    ##generation of list, handle DWD dataset
    station_list_dwd<-handle_dwd(action="list_stations",
                              location=c(9.40,48.49),
                              time_interval=c(19900101,20201231))

    Quality control summary produced by fix_weather()

    Station_name

    Region

    Station_ID

    Elevation

    Latitude

    Longitude

    Begin

    End

    Distance

    Lenningen-Schopfloch

    Baden-W<fc>rttemberg

    2949

    758

    48.5396

    9.5275

    19470801

    20020131

    10.92

    M<fc>nsingen-Apfelstetten

    Baden-W<fc>rttemberg

    3402

    755

    48.3851

    9.4837

    19460101

    20211231

    13.21

    Reutlingen-Betzingen

    Baden-W<fc>rttemberg

    4166

    360

    48.5051

    9.1833

    19880501

    20021231

    16.10

    N<fc>rtingen

    Baden-W<fc>rttemberg

    3671

    280

    48.6321

    9.3261

    19550101

    19910331

    16.72

    Trochtelfingen-Wilsingen

    Baden-W<fc>rttemberg

    5107

    769

    48.2814

    9.3137

    19880601

    20031231

    24.06

    Stuttgart-Echterdingen

    Baden-W<fc>rttemberg

    4931

    371

    48.6883

    9.2235

    19530101

    20211231

    25.61

  3. Identify suitable weather stations for patching gaps

    There are weather stations in the DWD list that are mostly complete I think, but for the sake of me trying to create a list (and we can never be sure that they indeed contain complete data with just the range of years shown), I will try to use 3 stations: 1, 3, and 6 from the weather station list.

  4. Download weather data for promising stations, convert them to chillR format and compile them in a list

    ##attempt 3
    require(reshape2)
    require(kableExtra)
    require(ggplot2)
    library(chillR)
    
    station_list_dwd<-handle_dwd(action="list_stations",location=c(9.40,48.49),time_interval=c(19900101,20201231))
    
    positions_in_station_list<-c(1,3,6)
    
    patch_weather<-list()
    for(i in 1:length(positions_in_station_list))
      {patch_weather[[i]]<-handle_dwd(handle_dwd(action="download_weather",
                                       location=station_list_dwd$Station_ID[positions_in_station_list[i]],time_interval=c(19900101,20211231)))[[1]]
    names(patch_weather)[i]<-station_list_dwd$Station_name[positions_in_station_list[i]]
    }
    
    save_temperature_scenarios(patch_weather,"data/", "patch_weather")
  5. Use the patch_daily_temperatures function to fill gaps

    Instead of Chapter 30, this exercise is the last bit I worked on because I had been wracking my head on how to carry this out when the patch_daily_temps seemed to only work perfectly fine with ‘gsod’ datafiles. So many hiccups because I could not also clean the ‘gsod’ raw files with handle_gsod for some reason. It always returns an empty file. That was the reason I had to use ‘dwd’ datafiles and move the location to Swabian Orchard paradise.

    After so many attempts, I was able to figure out that the ‘old and deprecated’ function patch_daily_temperatures would be able to handle ‘dwd’ datafiles well– at least in my case. Please remember that the patched data here is just for the sake of doing this exercise. I will be using the complete authentic data from dwd for my location for the next exercises. So, here’s the product of my many hours of sweat and tears:

    require(reshape2)
    require(kableExtra)
    require(ggplot2)
    library(chillR)
    ##recall previous list()
      station_list_dwd<-handle_dwd(action="list_stations",location=c(9.40,48.49),time_interval=c(19900101,20201231))
    
    positions_in_station_list<-c(1,3,6)
    
    patch_weather<-list()
    for(i in 1:length(positions_in_station_list))
      {patch_weather[[i]]<-handle_dwd(handle_dwd(action="download_weather",
                                       location=station_list_dwd$Station_ID[positions_in_station_list[i]],time_interval=c(19900101,20211231)))[[1]]
    names(patch_weather)[i]<-station_list_dwd$Station_name[positions_in_station_list[i]]
    }
    
    
    ##assign the dataset with gaps to object
    MApfelstetten<-read_tab("data/MApfelstetten_chillR_weather_gaps.csv")
    
    #first patching, I used the old function 'patch_daily_temperatures'
    
       patched<-patch_daily_temperatures(weather = MApfelstetten,
                                  patch_weather = patch_weather)
    
    ##to summarize statistics:
    #   > patched$statistics
    #$`Lenningen-Schopfloch`
    #     mean_bias stdev_bias filled gaps_remain
    #Tmin    -1.709      2.032     80        1027
    #Tmax    -0.085      1.283     80        1027
    
    #$`Reutlingen-Betzingen`
     #    mean_bias stdev_bias filled gaps_remain
    #Tmin    -2.062      1.664      0        1027
    #Tmax    -2.999      1.501      0        1027
    
    #$`Stuttgart-Echterdingen`
     #    mean_bias stdev_bias filled gaps_remain
    #Tmin    -1.950      1.595   1027           0
    #Tmax    -2.352      1.383   1027           0
    
    
    ##thus no need to limit mean bias to 1 and stdevbias to 2
    write.csv(patched$weather, "data/MApfelstetten_patched.csv", row.names = FALSE)
  6. Investigate the results - have all gaps been filled?

    library(chillR)
    patched<-read.csv("data/MApfelstetten_patched.csv")
    
    post_patch_stats<-fix_weather(patched)$QC
    
    write.csv(post_patch_stats, "data/MApfelstetten_post_patchstats.csv", row.names = FALSE)

    Data completeness table for the weather record from MApfelstetten, after applying the patch procedure

    Season

    End_year

    Season_days

    Data_days

    Missing_Tmin

    Missing_Tmax

    Incomplete_days

    Perc_complete

    1989/1990

    1990

    365

    365

    0

    0

    0

    100

    1990/1991

    1991

    365

    365

    0

    0

    0

    100

    1991/1992

    1992

    366

    366

    0

    0

    0

    100

    1992/1993

    1993

    365

    365

    0

    0

    0

    100

    1993/1994

    1994

    365

    365

    0

    0

    0

    100

    1994/1995

    1995

    365

    365

    0

    0

    0

    100

    1995/1996

    1996

    366

    366

    0

    0

    0

    100

    1996/1997

    1997

    365

    365

    0

    0

    0

    100

    1997/1998

    1998

    365

    365

    0

    0

    0

    100

    1998/1999

    1999

    365

    365

    0

    0

    0

    100

    1999/2000

    2000

    366

    366

    0

    0

    0

    100

    2000/2001

    2001

    365

    365

    0

    0

    0

    100

    2001/2002

    2002

    365

    365

    0

    0

    0

    100

    2002/2003

    2003

    365

    365

    0

    0

    0

    100

    2003/2004

    2004

    366

    366

    0

    0

    0

    100

    2004/2005

    2005

    365

    365

    0

    0

    0

    100

    2005/2006

    2006

    365

    365

    0

    0

    0

    100

    2006/2007

    2007

    365

    365

    0

    0

    0

    100

    2007/2008

    2008

    366

    366

    0

    0

    0

    100

    2008/2009

    2009

    365

    365

    0

    0

    0

    100

    2009/2010

    2010

    365

    365

    0

    0

    0

    100

    2010/2011

    2011

    365

    365

    0

    0

    0

    100

    2011/2012

    2012

    366

    366

    0

    0

    0

    100

    2012/2013

    2013

    365

    365

    0

    0

    0

    100

    2013/2014

    2014

    365

    365

    0

    0

    0

    100

    2014/2015

    2015

    365

    365

    0

    0

    0

    100

    2015/2016

    2016

    366

    366

    0

    0

    0

    100

    2016/2017

    2017

    365

    365

    0

    0

    0

    100

    2017/2018

    2018

    365

    365

    0

    0

    0

    100

    2018/2019

    2019

    365

    365

    0

    0

    0

    100

    2019/2020

    2020

    366

    366

    0

    0

    0

    100

  7. If necessary, repeat until you have a dataset you can work with in further analyses

    Since everything has been filled, no need to repeat procedure. What I wish to repeat though is that I will not use the produced dataset for further exercises.

11 Generating temperature scenarios

Exercise:

  1. For the location you chose for your earlier analyses, use chillR’s weather generator to produce 100 years of synthetic temperature data.

    I will be using Swabian Orchard Paradise (48.4931° N, 9.3998° E) for this portion.

    library(chillR)
    library(ggplot2)
    
    
    MApfelstetten_weather<-read.csv("data/MApfelstetten_weather.csv")
    
    Temp<-temperature_generation(MApfelstetten_weather,
                                 years=c(1998,2005),
                                 sim_years = c(2001,2100))
    
    Temperatures<-cbind(MApfelstetten_weather[
       which(MApfelstetten_weather$Year %in% 1998:2005),] ,Data_source="observed")
    
     Temperatures<-rbind(Temperatures,
                         cbind(Temp[[1]][,c("Year","Month","Day","Tmin","Tmax")],
                               Data_source="simulated"))
    
    Temperatures[,"Date"]<-as.Date(ISOdate(2000, Temperatures$Month, Temperatures$Day))
    
    
    # chill_comparison_full_seasons<-read.csv("data/chill_comparison_full_seasons.csv")
    
    ggplot(data=Temperatures, aes(Date,Tmin)) +
      geom_smooth(aes(colour = factor(Year))) +
      facet_wrap(vars(Data_source)) +
      theme_bw(base_size = 20) +
      theme(legend.position = "none") +
      scale_x_date(date_labels = "%b")
    
    ggplot(data=Temperatures, aes(Date,Tmax)) +
      geom_smooth(aes(colour = factor(Year))) +
      facet_wrap(vars(Data_source)) +
      theme_bw(base_size = 20) +
      theme(legend.position = "none") +
      scale_x_date(date_labels = "%b")

  2. Calculate winter chill (in Chill Portions) for your synthetic weather, and illustrate your results as histograms and cumulative distributions.

    library(chillR)
    library(ggplot2)
    
    MApfelstetten_weather<-read.csv("data/MApfelstetten_weather.csv")
    
    Temp<-temperature_generation(MApfelstetten_weather,
                                 years=c(1998,2005),
                                 sim_years = c(2001,2100))
    
    Temperatures<-cbind(MApfelstetten_weather[
       which(MApfelstetten_weather$Year %in% 1998:2005),] ,Data_source="observed")
    
     Temperatures<-rbind(Temperatures,
                         cbind(Temp[[1]][,c("Year","Month","Day","Tmin","Tmax")],
                               Data_source="simulated"))
    
    Temperatures[,"Date"]<-as.Date(ISOdate(2000, Temperatures$Month, Temperatures$Day))
    
    
     chill_observed<-chilling(
       stack_hourly_temps(
         Temperatures[which(Temperatures$Data_source=="observed"),],
         latitude = 48.49),
       Start_JDay = 305,
       End_JDay = 59)
    
     chill_simulated<-chilling(
       stack_hourly_temps(
         Temperatures[which(Temperatures$Data_source=="simulated"),],
         latitude = 48.49),
       Start_JDay = 305,
       End_JDay = 59)
    
     chill_comparison<-cbind(chill_observed ,Data_source="observed")
     chill_comparison<-rbind(chill_comparison,
                             cbind(chill_simulated ,Data_source="simulated"))
    
     chill_comparison_full_seasons<-chill_comparison[
       which(chill_comparison$Perc_complete==100),]
    
    ggplot(chill_comparison_full_seasons, aes(x=Chill_portions)) + 
      geom_histogram(binwidth=1,aes(fill = factor(Data_source))) +
      theme_bw(base_size = 20) +
      labs(fill = "Data source") +
      xlab("Chill accumulation (Chill Portions)") +
      ylab("Frequency")
    
    
    
    chill_simulations<-chill_comparison_full_seasons[
      which(chill_comparison_full_seasons$Data_source=="simulated"),]
    
    ggplot(chill_simulations, aes(x=Chill_portions)) +
      stat_ecdf(geom = "step",lwd=1.5,col="blue") +
      ylab("Cumulative probability") +
      xlab("Chill accumulation (in Chill Portions)") +
      theme_bw(base_size = 20)
    
    write.csv(chill_comparison_full_seasons,"data/chill_comparison_full_seasons.csv", row.names = FALSE)

    Season

    End_year

    Season_days

    Data_days

    Perc_complete

    Chilling_Hours

    Utah_Model

    Chill_portions

    GDH

    Data_source

    1998/1999

    1999

    120

    120

    100

    1016

    700.5

    57.31928

    497.6681

    observed

    1999/2000

    2000

    120

    120

    100

    1209

    810.0

    62.26994

    632.8378

    observed

    2000/2001

    2001

    121

    121

    100

    1621

    1262.5

    74.80941

    1027.6474

    observed

    2001/2002

    2002

    120

    120

    100

    1220

    992.5

    66.90289

    855.7386

    observed

    2002/2003

    2003

    120

    120

    100

    1301

    1055.0

    66.67006

    702.7373

    observed

    2003/2004

    2004

    120

    120

    100

    1211

    845.0

    65.66647

    1081.6018

    observed

    ##amount of chill that is exceeded in 90% of all years
    quantile(chill_simulations$Chill_portions, 0.1)
      #  10% 
    #59.0096 
    
    ##chill at 50% confidence interval (25th and 75th percentile)
    quantile(chill_simulations$Chill_portions, c(0.25,0.75))
       #  25%      75% 
    #62.81098 71.61968 
  3. Produce similar plots for the number of freezing hours (<0°C) in April (or October, if your site is in the Southern Hemisphere) for your location of interest.

    #attempt 3, freezing() function created based on chillR's chilling function
    library(chillR)
    library(ggplot2)
    
    ##synthetic weather generation 
    MApfelstetten_weather<-read.csv("data/MApfelstetten_weather.csv")
    
    Temp<-temperature_generation(MApfelstetten_weather,
                                 years=c(1998,2005),
                                 sim_years = c(2001,2100))
    
    Temperatures<-cbind(MApfelstetten_weather[
       which(MApfelstetten_weather$Year %in% 1998:2005),] ,Data_source="observed")
    
     Temperatures<-rbind(Temperatures,
                         cbind(Temp[[1]][,c("Year","Month","Day","Tmin","Tmax")],
                               Data_source="simulated"))
    
    Temperatures[,"Date"]<-as.Date(ISOdate(2000, Temperatures$Month, Temperatures$Day))
    
    ##creating function calculating freezing hours based on chillR function chilling()
    freezing<-function (hourtemps = NULL, Start_JDay = 1, End_JDay = 366, THourly = NULL, 
        misstolerance = 50) 
    {if (is.null(hourtemps) & !is.null(THourly)) 
            hourtemps <- THourly
        if ((length(names(hourtemps)) == 2) & ("hourtemps" %in% names(hourtemps))) {
            QC <- hourtemps$QC
            hourtemps <- hourtemps$hourtemps
        }
        else QC <- NULL
        if (Start_JDay < End_JDay) {
            hourtemps[which(hourtemps$JDay >= Start_JDay & hourtemps$JDay <= 
                End_JDay), "sea"] <- hourtemps[which(hourtemps$JDay >= 
                Start_JDay & hourtemps$JDay <= End_JDay), "Year"]
        }
        else {
            hourtemps[which(hourtemps$JDay >= Start_JDay), "sea"] <- hourtemps[which(hourtemps$JDay >= 
                Start_JDay), "Year"] + 1
            hourtemps[which(hourtemps$JDay <= End_JDay), "sea"] <- hourtemps[which(hourtemps$JDay <= 
                End_JDay), "Year"]
        }
        if (Start_JDay < End_JDay) {
            relevant_days <- Start_JDay:End_JDay
        }
        else {
            relevant_days <- c(Start_JDay:366, 1:End_JDay)
        }
        normal_lines <- which(!(hourtemps$JDay == Start_JDay & hourtemps$Hour == 
            0))
        normal_lines <- normal_lines[which(normal_lines > 1)]
        hourtemps <- hourtemps[which(!is.na(hourtemps[, "Temp"])), 
            ]
        ##range for freezing hour counting, creating columns
        FH_range <- which(hourtemps$Temp < 
            0)
         hourtemps[, "FH_weights"] <- 0
        hourtemps[FH_range, "FH_weights"] <- 1
        hourtemps[, "FH"] <- 0
    
      seasons <- as.numeric(unique(hourtemps$sea))
      seasons <- seasons[!is.na(seasons)]
    
    
        ##storing seasons
        chillout <- data.frame(Season = paste(seasons - 1, "/", seasons, 
            sep = ""), End_year = seasons)
        if (End_JDay >= Start_JDay) 
            dt <- sapply(seasons, function(x) difftime(ISOdate(x - 
                1, 12, 31) + End_JDay * 86400, ISOdate(x - 1, 12, 
                31) + Start_JDay * 86400))
        if (End_JDay < Start_JDay) 
            dt <- sapply(seasons, function(x) difftime(ISOdate(x - 
                1, 12, 31) + End_JDay * 86400, ISOdate(x - 2, 12, 
                31) + Start_JDay * 86400))
        chillout[, "Season_days"] <- dt + 1
        chillout[, "Data_days"] <- sapply(seasons, function(x) length(which(hourtemps$sea == 
            x))/24)
        for (sea in seasons) {
            if ("no_Tmin" %in% names(hourtemps) & "no_Tmax" %in% 
                names(hourtemps)) 
                chillout[which(chillout$End_year == sea), "Interpolated_days"] <- length(which(hourtemps[which(hourtemps$sea == 
                    sea), "no_Tmin"] | hourtemps[which(hourtemps$sea == 
                    sea), "no_Tmax"]))/24
            chillout[, "Perc_complete"] <- NA
    
            ##creation of column for summed freezing hours
            chillout[which(chillout$End_year == sea), "Freezing_Hours"] <- sum(hourtemps[which(hourtemps$sea == 
                sea), "FH_weights"])
        }
        if ("no_Tmin" %in% names(hourtemps) & "no_Tmax" %in% names(hourtemps)) 
            chillout[, "Perc_complete"] <- (chillout[, "Data_days"] - 
                chillout[, "Interpolated_days"])/chillout[, "Season_days"] * 
                100
        else chillout[, "Perc_complete"] <- chillout[, "Data_days"]/chillout[, 
            "Season_days"] * 100
        chillout <- chillout[which(chillout$Perc_complete >= 100 - 
            misstolerance), ]
        return(chillout)
    }
    
    
    freezing_observed<-freezing(
      stack_hourly_temps(
        Temperatures[which(Temperatures$Data_source=="observed"),],
        latitude = 48.49),
      Start_JDay = 91,
      End_JDay = 120)
    freezing_simulated<-freezing(
      stack_hourly_temps(
        Temperatures[which(Temperatures$Data_source=="simulated"),],
        latitude = 48.49),
      Start_JDay = 91,
      End_JDay = 120)
    
    freeze_comparison<-cbind(freezing_observed ,Data_source="observed")
    freeze_comparison<-rbind(freeze_comparison,
                            cbind(freezing_simulated ,Data_source="simulated"))
    
    freeze_comparison_full_seasons<-freeze_comparison[
      which(freeze_comparison$Perc_complete==100),]
    
    write.csv(freeze_comparison_full_seasons, "data/freeze_comparison_full_seasons.csv", row.names = FALSE)
    
    ##generating the Cumulative Freezing (sum of Freezing hours) vs Frequency plot
    ggplot(freeze_comparison_full_seasons, aes(x=Freezing_Hours)) + 
      geom_histogram(binwidth=1,aes(fill = factor(Data_source))) +
      theme_bw(base_size = 20) +
      labs(fill = "Data source") +
      xlab("Cumulative Freezing Hours") +
      ylab("Frequency")
    
    ##generating the Cumulative Freezing (sum of Freezing hours) vs Cumulative Probability plot
    
    freeze_simulations<-freeze_comparison_full_seasons[
      which(freeze_comparison_full_seasons$Data_source=="simulated"),]
    
    ggplot(freeze_simulations, aes(x=Freezing_Hours)) +
      stat_ecdf(geom = "step",lwd=1.5,col="blue") +
      ylab("Cumulative probability") +
      xlab("Cumulative Freezing Hours") +
      theme_bw(base_size = 20)

    Season

    End_year

    Season_days

    Data_days

    Perc_complete

    Freezing_Hours

    Data_source

    1997/1998

    1998

    30

    30

    100

    42

    observed

    1998/1999

    1999

    30

    30

    100

    48

    observed

    1999/2000

    2000

    30

    30

    100

    42

    observed

    2000/2001

    2001

    30

    30

    100

    117

    observed

    2001/2002

    2002

    30

    30

    100

    48

    observed

    2002/2003

    2003

    30

    30

    100

    142

    observed

    ##amount of chill that is exceeded in 90% of all years
    quantile(freeze_simulations$Freezing_Hours, 0.1)
    #10% 
    #7.9 
    
    ##chill at 50% confidence interval (25th and 75th percentile)
    quantile(freeze_simulations$Freezing_Hours, c(0.25,0.75))
    #25%   75% 
    #23.00 99.25 

13 Historic temperature scenarios

Exercise:

  1. For the location you chose for previous exercises, produce historic temperature scenarios representing several years of the historic record (your choice).

    After looking at the mean annual temperatures in Germany in a paper by Glaser et. al (2017), I saw that the temperatures only started to rise steeply during 1980s (same as what we had in the samples), so I decided to take a look at those same years too. Table from the paper shown below:

    require(chillR)
    require(ggplot2)
    require(reshape2)
    require(kableExtra)
    
    #download station list for time interval of interest
    station_list<-handle_dwd(action="list_stations",
                             location=c(9.40,48.49),
                             time_interval = c(19730101,20191231))
    
    station_list
    
    # download weather data from M<fc>nsingen-Apfelstetten (2nd in list) since it is still in the list and I've dealt with it before
    MApfelstetten_weather_raw<-handle_dwd(action="download_weather",
                                  location=station_list$Station_ID[2],
                                  time_interval = c(19730101,20191231),
                                  station_list = station_list)
    
    # convert weather data to chillR format
    MApfelstetten_weather<-handle_dwd(MApfelstetten_weather_raw)
    
    # check record for missing data
    MApfelstetten_historic<-fix_weather(MApfelstetten_weather[[1]])$QC
    
    write.csv(MApfelstetten_historic, "data/MApfelstetten_historic.csv", row.names = FALSE)
    
    ##there is no gap for the data set, proceeding to generating scenarios for 1980s
    MApfelstetten_temps<-MApfelstetten_weather[[1]]
    
    scenario_1980<-temperature_scenario_from_records(weather=MApfelstetten_temps,year=1980)
    
    scenario_1980$'1980'$data
    
    scenario_1980
    
    temps_1980<-temperature_generation(weather=MApfelstetten_temps, years=c(1973,2019),
                                       sim_years=c(2001,2100),
                                       temperature_scenario = scenario_1980)
    
    scenario_1996<-temperature_scenario_from_records(weather=MApfelstetten_temps,year=1996)
    
    relative_scenario<-temperature_scenario_baseline_adjustment(
      baseline=scenario_1996,
      temperature_scenario = scenario_1980)
    
    temps_1980<-temperature_generation(weather=MApfelstetten_temps, years=c(1973,2019),
                                       sim_years=c(2001,2100),
                                       temperature_scenario = relative_scenario)
    
    all_past_scenarios<-temperature_scenario_from_records(
      weather=MApfelstetten_temps,
      year=c(1980,1990,2000,2010))
    
    adjusted_scenarios<-temperature_scenario_baseline_adjustment(
      baseline=scenario_1996,
      temperature_scenario = all_past_scenarios)
    
    all_past_scenario_temps<-temperature_generation(
      weather=MApfelstetten_temps,
      years=c(1973,2019),
      sim_years=c(2001,2100),
      temperature_scenario = adjusted_scenarios)
    
    all_past_scenario_temps

    Season

    End_year

    Season_days

    Data_days

    Missing_Tmin

    Missing_Tmax

    Incomplete_days

    Perc_complete

    1972/1973

    1973

    365

    365

    0

    0

    0

    100

    1973/1974

    1974

    365

    365

    0

    0

    0

    100

    1974/1975

    1975

    365

    365

    0

    0

    0

    100

    1975/1976

    1976

    366

    366

    0

    0

    0

    100

    1976/1977

    1977

    365

    365

    0

    0

    0

    100

    1977/1978

    1978

    365

    365

    0

    0

    0

    100

  2. Produce chill distributions for these scenarios and plot them.

    require(chillR)
    require(ggplot2)
    require(reshape2)
    require(kableExtra)
    
    ##copying and pasting the code from previous step
    
    #download station list for time interval of interest
    station_list<-handle_dwd(action="list_stations",
                             location=c(9.40,48.49),
                             time_interval = c(19730101,20191231))
    
    station_list
    
    # download weather data from M<fc>nsingen-Apfelstetten (2nd in list) since it is still in the list and I've dealt with it before
    MApfelstetten_weather_raw<-handle_dwd(action="download_weather",
                                  location=station_list$Station_ID[2],
                                  time_interval = c(19730101,20191231),
                                  station_list = station_list)
    
    # convert weather data to chillR format
    MApfelstetten_weather<-handle_dwd(MApfelstetten_weather_raw)
    
    # check record for missing data
    MApfelstetten_historic<-fix_weather(MApfelstetten_weather[[1]])$QC
    
    write.csv(MApfelstetten_historic, "data/MApfelstetten_historic.csv", row.names = FALSE)
    
    ##there is no gap for the data set, proceeding to generating scenarios for 1980s
    MApfelstetten_temps<-MApfelstetten_weather[[1]]
    
    write.csv(MApfelstetten_temps, "data/MApfelstetten_temps.csv", row.names = FALSE)
    
    scenario_1980<-temperature_scenario_from_records(weather=MApfelstetten_temps,year=1980)
    
    scenario_1980$'1980'$data
    
    scenario_1980
    
    temps_1980<-temperature_generation(weather=MApfelstetten_temps, years=c(1973,2019),
                                       sim_years=c(2001,2100),
                                       temperature_scenario = scenario_1980)
    
    scenario_1996<-temperature_scenario_from_records(weather=MApfelstetten_temps,year=1996)
    
    relative_scenario<-temperature_scenario_baseline_adjustment(
      baseline=scenario_1996,
      temperature_scenario = scenario_1980)
    
    temps_1980<-temperature_generation(weather=MApfelstetten_temps, years=c(1973,2019),
                                       sim_years=c(2001,2100),
                                       temperature_scenario = relative_scenario)
    
    all_past_scenarios<-temperature_scenario_from_records(
      weather=MApfelstetten_temps,
      year=c(1980,1990,2000,2010))
    
    adjusted_scenarios<-temperature_scenario_baseline_adjustment(
      baseline=scenario_1996,
      temperature_scenario = all_past_scenarios)
    
    all_past_scenario_temps<-temperature_generation(
      weather=MApfelstetten_temps,
      years=c(1973,2019),
      sim_years=c(2001,2100),
      temperature_scenario = adjusted_scenarios)
    
    all_past_scenario_temps
    
    
    #chill accumulation calculation
    chill_hist_scenario_list<-tempResponse_daily_list(all_past_scenario_temps,latitude=48.49,Start_JDay = 305, End_JDay = 59)
    
    chill_hist_scenario_list
    
    #converting list to data.frame
    scenarios<-names(chill_hist_scenario_list)[1:4]
    
    all_scenarios<-chill_hist_scenario_list[[scenarios[1]]]
    all_scenarios[,"scenario"]<-as.numeric(scenarios[1])
    
    
    
    for (sc in scenarios[2:4])
     all_scenarios<-rbind(all_scenarios,
                          cbind(chill_hist_scenario_list[[sc]],scenario=as.numeric(sc)))
    
    all_scenarios<-all_scenarios[which(all_scenarios$Perc_complete==100),]
    
    #compute the actual 'observed' chill for comparison
    actual_chill<-tempResponse_daily_list(MApfelstetten_temps,latitude=48.49,
                            Start_JDay = 305,
                            End_JDay = 59)[[1]]
    actual_chill<-actual_chill[which(actual_chill$Perc_complete==100),]
    
    actual_chill
    
    ggplot(data=all_scenarios,aes(scenario,Chill_Portions,
                                  fill=factor(scenario))) +
      geom_violin() +
      ylab("Chill accumulation (Chill Portions)") +
      xlab("Scenario year") +
      theme_bw(base_size=12) +
      ylim(c(0,90)) +
      geom_point(data=actual_chill,
                 aes(End_year,Chill_Portions,fill="black"),
                 col="blue",show.legend = FALSE) +
      scale_fill_discrete(name="Scenario",
                          breaks = unique(all_scenarios$scenario)) 
    
    ggplot(data=all_scenarios,aes(scenario,GDH,
                                  fill=factor(scenario))) +
      geom_violin() +
      ylab("Heat accumulation (Growing Degree Hours)") +
      xlab("Scenario year") +
      theme_bw(base_size=12) +
      ylim(c(0,7000)) +
      geom_point(data=actual_chill,
                 aes(End_year,GDH,fill="black"),
                 col="blue",show.legend = FALSE) +
      scale_fill_discrete(name="Scenario",
                          breaks = unique(all_scenarios$scenario)) 
    
    
    ?runn_mean
    
    temperature_means<-data.frame(Year=min(MApfelstetten_temps$Year):max(MApfelstetten_temps$Year),
                                    Tmin=aggregate(MApfelstetten_temps$Tmin,FUN="mean",
                                                   by=list(MApfelstetten_temps$Year))[,2],
                                  Tmax=aggregate(MApfelstetten_temps$Tmax,FUN="mean",
                                                 by=list(MApfelstetten_temps$Year))[,2])
    temperature_means[,"runn_mean_Tmin"]<-runn_mean(temperature_means$Tmin,15)
    temperature_means[,"runn_mean_Tmax"]<-runn_mean(temperature_means$Tmax,15)
    
    Tmin_regression<-lm(Tmin~Year, temperature_means)
    temperature_means[,"regression_Tmin"]<-Tmin_regression$coefficients[1]+
      Tmin_regression$coefficients[2]*temperature_means$Year
    
    Tmax_regression<-lm(Tmax~Year, temperature_means)
    temperature_means[,"regression_Tmax"]<-Tmax_regression$coefficients[1]+
      Tmax_regression$coefficients[2]*temperature_means$Year
    
    
    ggplot(temperature_means,aes(Year, Tmin)) + 
      geom_point() + 
      geom_line(data=temperature_means,aes(Year, runn_mean_Tmin),lwd=2,col="blue") + 
      geom_line(data=temperature_means,aes(Year, regression_Tmin),lwd=2,col="red") +
      theme_bw(base_size=12) +
      ylab("Mean monthly minimum temperature (°C)")
    
    ggplot(temperature_means,aes(Year, Tmax)) + 
      geom_point() + 
      geom_line(data=temperature_means,aes(Year, runn_mean_Tmax),lwd=2,col="blue") + 
      geom_line(data=temperature_means,aes(Year, regression_Tmax),lwd=2,col="red") +
      theme_bw(base_size=12) +
      ylab("Mean monthly maximum temperature (°C)")

14 Future temperature scenarios

Exercise:

Analyze the historic and future impact of climate change on three agroclimatic metrics of your choice, for the location you’ve chosen for your earlier analyses.

install.packages("ggplot2")
require(chillR)
require(ggplot2)
require(reshape2)
require(kableExtra)


MApfelstetten_temps<-read_tab("data/MApfelstetten_temps.csv")


RCPs<-c("rcp45,rcp85")
Times<-c(2050,2090)


for(RCP in RCPs)
for(Time in Times)
    {start_year <- Time-15
      end_year <- Time+15
      clim_scen <-getClimateWizardData(
      c(longitude = 9.399,latitude = 48.493),
      RCP,
      start_year,
      end_year,
      temperature_generation_scenarios = TRUE,
      baseline =c(1975, 2005),
      metric = "monthly_min_max_temps",
      GCMs = "all")
    save_temperature_scenarios(clim_scen,
                               "data/ClimateWizard",
                               paste0("MApfelstetten_futures_",Time,"_",RCP))}


scenario_1990<-temperature_scenario_from_records(MApfelstetten_temps,1990)
scenario_1996<-temperature_scenario_from_records(MApfelstetten_temps,1996)
adjustment_scenario<-temperature_scenario_baseline_adjustment(scenario_1996,scenario_1990)


RCPs<-c("rcp85")
Times<-c(2090)

for(RCP in RCPs)
  for(Time in Times)
    {
     clim_scen<-load_ClimateWizard_scenarios(
       "data/climateWizard",
       paste0("MApfelstetten_futures_",2090,"_","rcp85"))
     clim_scen_adjusted<-
       temperature_scenario_baseline_adjustment(baseline_temperature_scenario = adjustment_scenario, temperature_scenario = clim_scen)
     Temps<-temperature_generation(
       weather=MApfelstetten_temps,
       years=c(1973,2019),
       sim_years=c(2001,2101),
       temperature_scenario = clim_scen_adjusted)

     save_temperature_scenarios(
       Temps,
       "data/Weather",
       paste0("MApfelstetten_",2090,"_","rcp85"))
  }

scenario_1990<-temperature_scenario_from_records(MApfelstetten_temps,1990)
scenario_1996<-temperature_scenario_from_records(MApfelstetten_temps,1996)
adjustment_scenario<-temperature_scenario_baseline_adjustment(scenario_1996,scenario_1990)

clim_scen<-load_ClimateWizard_scenarios(
       "data/climateWizard",
       paste0("MApfelstetten_futures_",2090,"_","rcp45"))

clim_scen
clim_scen_adjusted<-
      temperature_scenario_baseline_adjustment(baseline_temperature_scenario = adjustment_scenario, temperature_scenario = clim_scen)

Temps<-temperature_generation(
       weather=MApfelstetten_temps,
       years=c(1973,2019),
       sim_years=c(2001,2101),
       temperature_scenario = clim_scen_adjusted)

     save_temperature_scenarios(
       Temps,
       "data/Weather",
       paste0("MApfelstetten_",2090,"_","rcp45"))


RCPs<-c("rcp85")
Times<-c(2085)


for(RCP in RCPs)
  for(Time in Times)
    {
     clim_scen<-load_ClimateWizard_scenarios(
       "data/climateWizard",
       paste0("MApfelstetten_futures_",Time,"_",RCP))
     clim_scen_adjusted<-
       temperature_scenario_baseline_adjustment(
         baseline_temperature_scenario=adjustment_scenario,
         temperature_scenario=clim_scen)
     Temps<-temperature_generation(
       weather=MApfelstetten_temps,
       years=c(1973,2019),
       sim_years=c(2001,2101),
       temperature_scenario = clim_scen_adjusted)

     save_temperature_scenarios(
       Temps,
       "data/Weather",
       paste0("MApfelstetten_",Time,"_",RCP))
  }


##
all_past_scenarios<-temperature_scenario_from_records(
  weather=MApfelstetten_temps,
  year=c(1980,1990,2000,2010))

adjusted_scenarios<-temperature_scenario_baseline_adjustment(
  baseline=scenario_1996,
  temperature_scenario = all_past_scenarios)

all_past_scenario_temps<-temperature_generation(
  weather=MApfelstetten_temps,
  years=c(1973,2019),
  sim_years=c(2001,2101),
  temperature_scenario = adjusted_scenarios)

save_temperature_scenarios(
  all_past_scenario_temps,
  "data/Weather",
  "MApfelstetten_historic")

frost_model<-function(x) step_model(x,data.frame(
  lower=c(-1000,0),
  upper=c(0,1000),
  weight=c(1,0)))
  
  models<-list(Chill_CP=Dynamic_Model,
             Heat_GDH=GDH,
             Frost_H=frost_model)
  
  Temps<-load_temperature_scenarios("data/Weather","MApfelstetten_historic")
chill_past_scenarios<-tempResponse_daily_list(
  Temps,
  latitude=48.49,
  Start_JDay = 305,
  End_JDay = 59,
  models=models,
  misstolerance = 10)
chill_observed<-tempResponse_daily_list(
  MApfelstetten_temps,
  latitude=48.49,
  Start_JDay = 305,
  End_JDay = 59,
  models=models,
  misstolerance = 10)
save_temperature_scenarios(chill_past_scenarios,
                           "data/chill",
                           "MApfelstetten_historic")
save_temperature_scenarios(chill_observed,
                           "data/chill",
                           "MApfelstetten_observed")
chill_past_scenarios<-load_temperature_scenarios(
  "data/chill",
  "MApfelstetten_historic")
chill_observed<-load_temperature_scenarios(
  "data/chill",
  "MApfelstetten_observed")

chills <-make_climate_scenario(
  chill_past_scenarios,
  caption = "Historic",
  historic_data = chill_observed,
  time_series = TRUE)

plot_climate_scenarios(
  climate_scenario_list=chills,
  metric="Chill_CP",
  metric_label="Chill (Chill Portions)")


RCPs<-c("rcp45", "rcp85")
Times<-c(2050,2090)

for(RCP in RCPs)
  for(Time in Times)
    {
    Temps<-load_temperature_scenarios(
      "data/Weather",
      paste0("MApfelstetten_",Time,"_",RCP))
    chill<-tempResponse_daily_list(
      Temps,
      latitude=48.49,
      Start_JDay = 305,
      End_JDay = 59,
      models=models,
      misstolerance = 10)
    save_temperature_scenarios(
      chill,
      "data/chill",
      paste0("MApfelstetten_",Time,"_",RCP))
  }


for(RCP in RCPs)
  for(Time in Times)
    {
    chill<-load_temperature_scenarios(
      "data/chill",
      paste0("MApfelstetten_",Time,"_",RCP))
    if(RCP=="rcp45") RCPcaption <- "RCP4.5"
    if(RCP=="rcp85") RCPcaption <- "RCP8.5"
    if(Time=="2050") Time_caption <- "2050"
    if(Time=="2090") Time_caption <- "2090"
    chills <-make_climate_scenario(
      chill,
      caption =c(RCPcaption, Time_caption),
      add_to = chills)
  }

#plotting
plot_climate_scenarios(
  climate_scenario_list=chills,
  metric="Chill_CP",
  metric_label="Chill (Chill Portions)",
  texcex=1.5)

plot_climate_scenarios(
    climate_scenario_list=chills,
   metric="Frost_H",
    metric_label="Frost hours",
   texcex=1.5)

Chill Portions- the Swabian Orchard Paradise will not have to worry about not being able to grow their trees up to 2090, it seems, even with RCP8.5 projection.

Heat (GDH)- per the projection, it will become warmer, thus can give opportunity for farmers to venture into more diverse tree species that require higher GDH in the future.

Frost hours- hours of frost is projected to decline, which is also good because less frost damage to flowers during critical times.

In general, these plots come in agreement with what we are expecting: that the location (which is in Germany) will be more fit for production- this was in line with what lecturers from other courses were saying. Global North will be benefiting more from climate change than the Global South.

15 Plotting future scenarios

Exercise: Produce similar plots for the weather station you selected for earlier exercises.

##install.packages("ggpmisc")
library(chillR)
library(ggplot2)
library(reshape2)
library(kableExtra)
library(ggpmisc)
library(patchwork)



chill_past_scenarios<-load_temperature_scenarios(
  "data/chill",
  "MApfelstetten_historic")
chill_observed<-load_temperature_scenarios(
  "data/chill",
  "MApfelstetten_observed")

chills <-make_climate_scenario(
  chill_past_scenarios,
  caption = "Historic",
  historic_data = chill_observed,
  time_series = TRUE)

RCPs<-c("rcp45","rcp85")
Times<-c(2050,2090)

for(RCP in RCPs)
  for(Time in Times)
  { chill<-load_temperature_scenarios(
    "data/chill",
    paste0("MApfelstetten_",Time,"_",RCP))
    if(RCP=="rcp45") RCPcaption <- "RCP4.5"
    if(RCP=="rcp85") RCPcaption <- "RCP8.5"
    if(Time=="2050") Time_caption <- "2050"
    if(Time=="2090") Time_caption <- "2090"
    chills <-make_climate_scenario(
      chill,
      caption =c(RCPcaption, Time_caption),
      add_to = chills)
  }

plot_climate_scenarios(
  climate_scenario_list=chills,
  metric="Frost_H",
  metric_label="Frost hours",
  texcex=1.5)

# We'll first process the past scenarios (element 1 of the chills list).
# Within the data element, we have a list of multiple data.frames for
# the various past scenarios.
# Using a 'for' loop, we cycle through all these data.frames.

for(nam in names(chills[[1]]$data))
  {
   # Extract the data frame.
   ch<-chills[[1]]$data[[nam]]
   # Add columns for the new information we have to add and fill them.
   ch[,"GCM"]<-"none"
   ch[,"RCP"]<-"none"
   ch[,"Year"]<-as.numeric(nam)
   
   # Now check if this is the first time we've gone through this loop.
   # If this is the first time, the ch data.frame becomes the output
   # object (past_simulated).
   # If it is not the first time ('else'), we add the current data.frame
   # to the 'past_simulated' object
  if(nam==names(chills[[1]]$data)[1])
    past_simulated<-ch else
      past_simulated<-rbind(past_simulated,ch)
  }

# We add another column called 'Scenario' and label all rows as 'Historic' 
past_simulated["Scenario"] <- "Historic"

kable(past_simulated[1:5,])  %>%
  kable_styling("striped", position = "left",font_size = 8)

# We'll want to add the historic observation too, so let's simplify the
# pointer to this information for easier use later

past_observed <- chills[[1]][["historic_data"]]

kable(past_observed[1:5,])  %>%
  kable_styling("striped", position = "left",font_size = 8)


# Extract future data
for(i in 2:length(chills))
  for(nam in names(chills[[i]]$data))
    {ch<-chills[[i]]$data[[nam]]
     ch[,"GCM"]<-nam
     ch[,"RCP"]<-chills[[i]]$caption[1]
     ch[,"Year"]<-chills[[i]]$caption[2]
     if(i==2&nam==names(chills[[i]]$data)[1])
       future_data<-ch else
         future_data<-rbind(future_data,ch)
  }

kable(future_data[1:5,])  %>%
  kable_styling("striped", position = "left",font_size = 8)


# Extract the model names
Models <- unique(future_data$GCM)
  
metric<-"Heat_GDH"
axis_label<-"Heat (in GDH)"

# get extreme values for the axis scale

rng<-range(past_observed[[metric]],
           past_simulated[[metric]],
           future_data[[metric]])  
rng


past_plot<-ggplot() +
  geom_boxplot(data = past_simulated,
               aes_string("as.numeric(Year)",metric,group="Year"),
               fill="skyblue")

past_plot



past_plot<-past_plot +
  scale_y_continuous(limits = c(0, round(rng[2] + rng[2]/10))) +
  labs(x = "Year", y = axis_label)

past_plot


past_plot<-past_plot +
  facet_grid(~ Scenario) +
  theme_bw(base_size = 15) 
  
past_plot

  
past_plot<-past_plot +  
  theme(strip.background = element_blank(),
        strip.text = element_text(face = "bold"),
        axis.text.x = element_text(angle=45, hjust=1)) 

past_plot



# add historic data
past_plot <- past_plot +
  geom_point(data = past_observed,
             aes_string("End_year",metric),
             col="blue")

past_plot



y<-2050

future_2050<-ggplot(data= future_data[which(future_data$Year==y),]) +
      geom_boxplot(aes_string("GCM", metric, fill="GCM"))

future_2050



future_2050<-future_2050 +
  facet_wrap(vars(RCP)) +
   scale_x_discrete(labels = NULL, expand = expansion(add = 1)) 


##stopped here
future_2050<-future_2050 +
  scale_y_continuous(limits = c(0, round(1.1*rng[2]))) +
    geom_text_npc(aes(npcx = "center", npcy = "top", label = Year), size = 5)

future_2050


future_2050<-future_2050 +
  theme_bw(base_size = 15) +
  theme(axis.ticks.y = element_blank(),
        axis.text = element_blank(),
        axis.title = element_blank(),
                 legend.position = "bottom",
                 legend.margin = margin(0, 0, 0, 0, "cm"),
                 legend.background = element_rect(),
                 strip.background = element_blank(),
                 strip.text = element_text(face = "bold"),
                 legend.box.spacing = unit(0, "cm"),
                 plot.subtitle = element_text(hjust = 0.5, vjust = -1, size = 15 * 1.05,
                                                       face = "bold")) 

future_2050


future_plot_list<-list()

for(y in c(2050,2090))
{
  future_plot_list[[which(y == c(2050,2090))]] <-
    ggplot(data= future_data[which(future_data$Year==y),]) +
    geom_boxplot(aes_string("GCM", metric, fill="GCM")) +
    facet_wrap(vars(RCP)) +
    scale_x_discrete(labels = NULL, expand = expansion(add = 1)) +
    scale_y_continuous(limits = c(0, round(round(1.1*rng[2])))) +
    geom_text_npc(aes(npcx = "center", npcy = "top", label = Year),
                  size = 5) +
    theme_bw(base_size = 15) +
    theme(axis.ticks.y = element_blank(),
          axis.text = element_blank(),
          axis.title = element_blank(),
          legend.position = "bottom",
          legend.margin = margin(0, 0, 0, 0, "cm"),
          legend.background = element_rect(),
          strip.background = element_blank(),
          strip.text = element_text(face = "bold"),
          legend.box.spacing = unit(0, "cm"),
          plot.subtitle = element_text(
            hjust = 0.5,
            vjust = -1,
            size = 15 * 1.05,
            face = "bold")) 
}

future_plot_list

library(patchwork)

both_plots<-past_plot + future_plot_list

both_plots



plot <- both_plots +
           plot_layout(guides = "collect",
                       widths = c(1,rep(1.8,length(future_plot_list))))
         


plot<-plot & theme(legend.position = "bottom",
                 legend.text = element_text(size=8),
                 legend.title = element_text(size=10),
                 axis.title.x=element_blank())


plot

         
metric<-"Chill_CP"
axis_label<-"Chill (in CP)"

# get extreme values for the axis scale

rng<-range(past_observed[[metric]],past_simulated[[metric]],future_data[[metric]])  
past_plot<-ggplot() +
  geom_boxplot(data = past_simulated, aes_string("as.numeric(Year)",metric,group="Year"),fill="skyblue") +
  scale_y_continuous(limits = c(0, round(round(1.1*rng[2])))) +
  labs(x = "Year", y = axis_label) +
  facet_grid(~ Scenario) +
  theme_bw(base_size = 15) +  
  theme(strip.background = element_blank(),
           strip.text = element_text(face = "bold"),
           axis.text.x = element_text(angle=45, hjust=1)) +
  geom_point(data = past_observed,aes_string("End_year",metric),col="blue")

future_plot_list<-list()

for(y in c(2050,2090))
{
  future_plot_list[[which(y == c(2050,2090))]] <- ggplot(data= future_data[which(future_data$Year==y),]) +
      geom_boxplot(aes_string("GCM", metric, fill="GCM")) +
  facet_wrap(vars(RCP)) +
   scale_x_discrete(labels = NULL, expand = expansion(add = 1)) +
  scale_y_continuous(limits = c(0, round(round(1.1*rng[2])))) +
    geom_text_npc(aes(npcx = "center", npcy = "top", label = Year), size = 5) +
  theme_bw(base_size = 15) +
  theme(axis.ticks.y = element_blank(),
        axis.text = element_blank(),
        axis.title = element_blank(),
                 legend.position = "bottom",
                 legend.margin = margin(0, 0, 0, 0, "cm"),
                 legend.background = element_rect(),
                 strip.background = element_blank(),
                 strip.text = element_text(face = "bold"),
                 legend.box.spacing = unit(0, "cm"),
                 plot.subtitle = element_text(hjust = 0.5, vjust = -1, size = 15 * 1.05,
                                                       face = "bold")) 
}

plot<- (past_plot +
        future_plot_list +
        plot_layout(guides = "collect",
                       widths = c(1,rep(1.8,length(future_plot_list))))
       ) & theme(legend.position = "bottom",
                 legend.text = element_text(size=8),
                 legend.title = element_text(size=10),
                 axis.title.x=element_blank())
plot

         
metric<-"Frost_H"
axis_label<-"Frost duration (in hours)"

# get extreme values for the axis scale

rng<-range(past_observed[[metric]],past_simulated[[metric]],future_data[[metric]])  
past_plot<-ggplot() +
  geom_boxplot(data = past_simulated, aes_string("as.numeric(Year)",metric,group="Year"),fill="skyblue") +
  scale_y_continuous(limits = c(0, round(round(1.1*rng[2])))) +
  labs(x = "Year", y = axis_label) +
  facet_grid(~ Scenario) +
  theme_bw(base_size = 15) +  
  theme(strip.background = element_blank(),
           strip.text = element_text(face = "bold"),
           axis.text.x = element_text(angle=45, hjust=1)) +
  geom_point(data = past_observed,aes_string("End_year",metric),col="blue")

future_plot_list<-list()

for(y in c(2050,2090))
{
  future_plot_list[[which(y == c(2050,2090))]] <- ggplot(data= future_data[which(future_data$Year==y),]) +
      geom_boxplot(aes_string("GCM", metric, fill="GCM")) +
  facet_wrap(vars(RCP)) +
   scale_x_discrete(labels = NULL, expand = expansion(add = 1)) +
  scale_y_continuous(limits = c(0, round(round(1.1*rng[2])))) +
    geom_text_npc(aes(npcx = "center", npcy = "top", label = Year), size = 5) +
  theme_bw(base_size = 15) +
  theme(axis.ticks.y = element_blank(),
        axis.text = element_blank(),
        axis.title = element_blank(),
                 legend.position = "bottom",
                 legend.margin = margin(0, 0, 0, 0, "cm"),
                 legend.background = element_rect(),
                 strip.background = element_blank(),
                 strip.text = element_text(face = "bold"),
                 legend.box.spacing = unit(0, "cm"),
                 plot.subtitle = element_text(hjust = 0.5, vjust = -1, size = 15 * 1.05,
                                                       face = "bold")) 
}

plot<- (past_plot +
        future_plot_list +
        plot_layout(guides = "collect",
                       widths = c(1,rep(1.8,length(future_plot_list))))
       ) & theme(legend.position = "bottom",
                 legend.text = element_text(size=8),
                 legend.title = element_text(size=10),
                 axis.title.x=element_blank())
plot



plot_scenarios_gg<-function(past_observed,
                            past_simulated,
                            future_data,
                            metric,
                            axis_label)
{
  rng<-range(past_observed[[metric]],
             past_simulated[[metric]],
             future_data[[metric]])  
  past_plot<-ggplot() +
    geom_boxplot(data = past_simulated,
                 aes_string("as.numeric(Year)",metric,group="Year"),
                 fill="skyblue") +
    scale_y_continuous(limits = c(0, round(round(1.1*rng[2])))) +
    labs(x = "Year", y = axis_label) +
    facet_grid(~ Scenario) +
    theme_bw(base_size = 15) +  
    theme(strip.background = element_blank(),
          strip.text = element_text(face = "bold"),
          axis.text.x = element_text(angle=45, hjust=1)) +
    geom_point(data = past_observed,
               aes_string("End_year",metric),
               col="blue")
  
  future_plot_list<-list()
  
  for(y in c(2050,2090))
  {
    future_plot_list[[which(y == c(2050,2090))]] <-
      ggplot(data= future_data[which(future_data$Year==y),]) +
      geom_boxplot(aes_string("GCM", metric, fill="GCM")) +
      facet_wrap(vars(RCP)) +
      scale_x_discrete(labels = NULL, expand = expansion(add = 1)) +
      scale_y_continuous(limits = c(0, round(round(1.1*rng[2])))) +
      geom_text_npc(aes(npcx = "center", npcy = "top", label = Year),
                    size = 5) +
      theme_bw(base_size = 15) +
      theme(axis.ticks.y = element_blank(),
            axis.text = element_blank(),
            axis.title = element_blank(),
            legend.position = "bottom",
            legend.margin = margin(0, 0, 0, 0, "cm"),
            legend.background = element_rect(),
            strip.background = element_blank(),
            strip.text = element_text(face = "bold"),
            legend.box.spacing = unit(0, "cm"),
            plot.subtitle = element_text(hjust = 0.5,
                                         vjust = -1,
                                         size = 15 * 1.05,
                                         face = "bold")) 
  }
  
  plot<- (past_plot +
            future_plot_list +
            plot_layout(guides = "collect",
                        widths = c(1,rep(1.8,length(future_plot_list))))
        ) & theme(legend.position = "bottom",
                  legend.text = element_text(size=8),
                  legend.title = element_text(size=10),
                  axis.title.x=element_blank())
  plot
  
}





##plotting and saving plots (I manually created the "plots" folder)
plot_scenarios_gg(past_observed=past_observed,
                  past_simulated=past_simulated,
                  future_data=future_data,
                  metric="Heat_GDH",
                  axis_label="Heat (in Growing Degree Hours)")
ggsave("data/plots/Heat_GDH.png", width = 20,
       height = 10, units = "cm",
       dpi = 600)

plot_scenarios_gg(past_observed=past_observed,
                  past_simulated=past_simulated,
                  future_data=future_data,
                  metric="Chill_CP",
                  axis_label="Chill (in Chill Portions)")

ggsave("data/plots/Chill_CP.png", width = 20,
       height = 10, units = "cm",
       dpi = 600)

plot_scenarios_gg(past_observed=past_observed,
                  past_simulated=past_simulated,
                  future_data=future_data,
                  metric="Frost_H",
                  axis_label="Frost duration (in hours)")

ggsave("data/plots/Frost_H.png", width = 20,
       height = 10, units = "cm",
       dpi = 600)

16 Chill model comparison

Exercise:

  1. Perform a similar analysis for the location you’ve chosen for your exercises.

    I just want to share that I got super stuck in installing dormancyR because typing in the devtools library constantly gave an error that looked like this:

    ##Error: package or namespace load failed for ‘devtools’ in loadNamespace(i, c(lib.loc, .libPaths()), versionCheck = vI[[i]]):
    ## namespace ‘htmltools’ 0.5.3 is already loaded, but >= 0.5.4 is required

    I had to install and use the “remotes” package for it to work. Anyway, here’s the code:

    library(chillR)
    ##install.packages("remotes")
    ##install.packages("colorRamps")
    install.packages("gganimate")
    library(remotes)
    remotes::install_github("EduardoFernandezC/dormancyR")
    library(dormancyR)
    library(ggplot2)
    library(reshape2)
    library(colorRamps)
    library(kableExtra)
    library(patchwork)
    library(gganimate)
    
    RCPs<-c("rcp45","rcp85")
    Times<-c(2050,2090)
    
    hourly_models <- list(Chilling_units = chilling_units,
         Low_chill = low_chill_model,
         Modified_Utah = modified_utah_model,
         North_Carolina = north_carolina_model,
         Positive_Utah = positive_utah_model,
         Chilling_Hours = Chilling_Hours,
         Utah_Chill_Units = Utah_Model,
         Chill_Portions = Dynamic_Model)
    daily_models<-list(Rate_of_Chill = rate_of_chill,
                       Chill_Days = chill_days,
                       Exponential_Chill = exponential_chill,
                       Triangula_Chill_Haninnen = triangular_chill_1,
                       Triangular_Chill_Legave = triangular_chill_2)
    
    metrics<-c(names(daily_models),names(hourly_models))
    
    model_labels=c("Rate of Chill",
                   "Chill Days",
                   "Exponential Chill",
                   "Triangular Chill (Häninnen)",
                   "Triangular Chill (Legave)",
                   "Chilling Units",
                   "Low-Chill Chill Units",
                   "Modified Utah Chill Units",
                   "North Carolina Chill Units",
                   "Positive Utah Chill Units",
                   "Chilling Hours",
                   "Utah Chill Units",
                   "Chill Portions")
    
    kable(data.frame(Metric=model_labels,'Function name'=metrics))  %>%
      kable_styling("striped", position = "left",font_size = 10)
    
    MApfelstetten_temps<-read_tab("data/MApfelstetten_temps.csv")
    Temps<-load_temperature_scenarios("data/Weather","MApfelstetten_historic")
    
    
    Start_JDay<-305
    End_JDay<-59
    
    daily_models_past_scenarios<-tempResponse_list_daily(
      Temps,
      Start_JDay = Start_JDay,
      End_JDay = End_JDay,
      models=daily_models)
    
    
    hourly_models_past_scenarios<-tempResponse_daily_list(
      Temps,
      latitude=48.493,
      Start_JDay = Start_JDay,
      End_JDay = End_JDay,
      models=hourly_models,
      misstolerance = 10)
    
    #past_scenarios<-daily_models_past_scenarios
    past_scenarios<-lapply(
      names(daily_models_past_scenarios),
      function(x)
        cbind(daily_models_past_scenarios[[x]],
              hourly_models_past_scenarios[[x]][,names(hourly_models)]))
    names(past_scenarios)<-names(daily_models_past_scenarios)
    
    daily_models_observed<-tempResponse_daily(
      MApfelstetten_temps,
      Start_JDay = Start_JDay,
      End_JDay = End_JDay,
      models=daily_models)
    daily_models_observed<-
      daily_models_observed[which(daily_models_observed$Perc_complete>90),]
    hourly_models_observed<-tempResponse_daily_list(
      MApfelstetten_temps,
      latitude=48.493,
      Start_JDay = Start_JDay,
      End_JDay = End_JDay,
      models=hourly_models,
      misstolerance = 10)
    
    past_observed<-cbind(
      daily_models_observed,
      hourly_models_observed[[1]][,names(hourly_models)])
    
    save_temperature_scenarios(past_scenarios,
                               "data/chill",
                               "MApfelstetten_multichill_historic")
    write.csv(past_observed,
              "data/chill/MApfelstetten_multichill_observed.csv",
              row.names=FALSE)
    
    
    ##same procedure for future data
    RCPs<-c("rcp45","rcp85")
    Times<-c(2050,2090)
    
    for(RCP in RCPs)
      for(Time in Times)
        {
        Temps<-load_temperature_scenarios(
          "data/Weather",
          paste0("MApfelstetten_",Time,"_",RCP))
    
        daily_models_future_scenarios<-tempResponse_list_daily(
          Temps,
          Start_JDay = Start_JDay,
          End_JDay = End_JDay,
          models=daily_models)
        daily_models_future_scenarios<-lapply(
          daily_models_future_scenarios,
          function(x) x[which(x$Perc_complete>90),])
    
        hourly_models_future_scenarios<-
          tempResponse_daily_list(
            Temps,
            latitude=48.493,
            Start_JDay = Start_JDay,
            End_JDay = End_JDay,
            models=hourly_models,
            misstolerance = 10)
    
        future_scenarios<-daily_models_future_scenarios
        future_scenarios<-lapply(
          names(future_scenarios),
          function(x)
            cbind(future_scenarios[[x]],
                  hourly_models_future_scenarios[[x]][,names(hourly_models)]))
        names(future_scenarios)<-names(daily_models_future_scenarios)
    
        chill<-future_scenarios
        save_temperature_scenarios(
          chill,
          "data/chill",
          paste0("MApfelstetten_multichill_",Time,"_",RCP))
    }
    
    
    ##loading and producing chill scenarios
    chill_past_scenarios<-load_temperature_scenarios(
      "data/chill",
      "MApfelstetten_multichill_historic")
    chill_observed<-read_tab("data/chill/MApfelstetten_multichill_observed.csv")
    
    chills <-make_climate_scenario(chill_past_scenarios,
                                   caption = "Historic",
                                   historic_data = chill_observed,
                                   time_series = TRUE)
    
    for(RCP in RCPs)
      for(Time in Times)
        {
        chill<-load_temperature_scenarios(
          "data/chill",
          paste0("MApfelstetten_multichill_",Time,"_",RCP))
        if(RCP=="rcp45") RCPcaption <- "RCP4.5"
        if(RCP=="rcp85") RCPcaption <- "RCP8.5"
        if(Time=="2050") Time_caption <- "2050"
        if(Time=="2090") Time_caption <- "2090"
        chills <-make_climate_scenario(chill,
                                       caption =c(RCPcaption,Time_caption),
                                       add_to = chills)
    }
    
    
    ##heat map for Safe Winter chill
    for(i in 1:length(chills))
       {ch<-chills[[i]]
        if(ch$caption[1]=="Historic")
          {GCMs<-rep("none",length(names(ch$data)))
           RCPs<-rep("none",length(names(ch$data)))
           Years<-as.numeric(ch$labels)
           Scenario<-rep("Historic",length(names(ch$data)))} else
             {GCMs<-names(ch$data)
              RCPs<-rep(ch$caption[1],length(names(ch$data)))
              Years<-rep(as.numeric(ch$caption[2]),length(names(ch$data)))
              Scenario<-rep("Future",length(names(ch$data)))}
        for(nam in names(ch$data))
        {for(met in metrics)
          {temp_res<-data.frame(Metric=met,
                       GCM=GCMs[which(nam==names(ch$data))],
                       RCP=RCPs[which(nam==names(ch$data))],
                       Year=Years[which(nam==names(ch$data))],
                       Result=quantile(ch$data[[nam]][,met],0.1),
                       Scenario=Scenario[which(nam==names(ch$data))])
          if(i==1&nam==names(ch$data)[1]&met==metrics[1])
            results<-temp_res else
              results<-rbind(results,temp_res)
    
        }}}
    
    for(met in metrics)
      results[which(results$Metric==met),"SWC"]<-
        results[which(results$Metric==met),"Result"]/
          results[which(results$Metric==met&results$Year==1980),"Result"]-1
    
    
    
    rng = range(results$SWC)
    
    p_future<-ggplot(results[which(!results$GCM=="none"),],
                     aes(GCM, y=factor(Metric, levels=metrics),
                         fill = SWC)) +
      geom_tile()
    
    p_future
    
    
    p_future <-
      p_future +
      facet_grid(RCP ~ Year) 
    p_future
    
    #make text smaller
    p_future <-
      p_future +
      theme_bw(base_size = 15) +
      theme(axis.text = element_text(size=8))
    
    p_future
    
    library(colorRamps)
    p_future <-
      p_future +
      scale_fill_gradientn(colours=matlab.like(15),
                           labels = scales::percent,
                           limits=rng)
    p_future
    
    
    p_future <-
      p_future  +
      theme(axis.text.x = element_text(angle = 75, hjust = 1, vjust = 1)) +
      labs(fill = "Change in\nSafe Winter Chill\nsince 1975") +
      scale_y_discrete(labels=model_labels) +
      ylab("Chill metric")
    p_future
    
    
    #plot of past scenarios
    
    p_past<-
      ggplot(results[which(results$GCM=="none"),],
             aes(Year, y=factor(Metric, levels=metrics),
                 fill = SWC)) +
      geom_tile()
    
    p_past<-
      p_past +
      theme_bw(base_size = 12) +
      theme(axis.text = element_text(size=8))
    
    p_past<-
      p_past +
      scale_fill_gradientn(colours=matlab.like(15),
                           labels = scales::percent,
                           limits=rng)
    
    p_past<-
      p_past +
      scale_x_continuous(position = "top") 
    
    p_past<-
      p_past +
      labs(fill = "Change in\nSafe Winter Chill\nsince 1975") +
      scale_y_discrete(labels=model_labels) +
      ylab("Chill metric")
    
    
    chill_comp_plot<-
      (p_past +
         p_future +
         plot_layout(guides = "collect",nrow=2, heights=c(1,2))) &
      theme(legend.position = "right",strip.background = element_blank(),
            strip.text = element_text(face = "bold"))
    
    chill_comp_plot
  2. Make a heat map illustrating past and future changes in Safe Winter Chill, relative to a past scenario, for the 13 chill models used here

    library(chillR)
    ##install.packages("remotes")
    ##install.packages("colorRamps")
    install.packages("gganimate")
    library(remotes)
    remotes::install_github("EduardoFernandezC/dormancyR")
    library(dormancyR)
    library(ggplot2)
    library(reshape2)
    library(colorRamps)
    library(kableExtra)
    library(patchwork)
    library(gganimate)
    
    RCPs<-c("rcp45","rcp85")
    Times<-c(2050,2090)
    
    hourly_models <- list(Chilling_units = chilling_units,
         Low_chill = low_chill_model,
         Modified_Utah = modified_utah_model,
         North_Carolina = north_carolina_model,
         Positive_Utah = positive_utah_model,
         Chilling_Hours = Chilling_Hours,
         Utah_Chill_Units = Utah_Model,
         Chill_Portions = Dynamic_Model)
    daily_models<-list(Rate_of_Chill = rate_of_chill,
                       Chill_Days = chill_days,
                       Exponential_Chill = exponential_chill,
                       Triangula_Chill_Haninnen = triangular_chill_1,
                       Triangular_Chill_Legave = triangular_chill_2)
    
    metrics<-c(names(daily_models),names(hourly_models))
    
    model_labels=c("Rate of Chill",
                   "Chill Days",
                   "Exponential Chill",
                   "Triangular Chill (Häninnen)",
                   "Triangular Chill (Legave)",
                   "Chilling Units",
                   "Low-Chill Chill Units",
                   "Modified Utah Chill Units",
                   "North Carolina Chill Units",
                   "Positive Utah Chill Units",
                   "Chilling Hours",
                   "Utah Chill Units",
                   "Chill Portions")
    
    kable(data.frame(Metric=model_labels,'Function name'=metrics))  %>%
      kable_styling("striped", position = "left",font_size = 10)
    
    MApfelstetten_temps<-read_tab("data/MApfelstetten_temps.csv")
    Temps<-load_temperature_scenarios("data/Weather","MApfelstetten_historic")
    
    
    Start_JDay<-305
    End_JDay<-59
    
    daily_models_past_scenarios<-tempResponse_list_daily(
      Temps,
      Start_JDay = Start_JDay,
      End_JDay = End_JDay,
      models=daily_models)
    
    
    hourly_models_past_scenarios<-tempResponse_daily_list(
      Temps,
      latitude=48.493,
      Start_JDay = Start_JDay,
      End_JDay = End_JDay,
      models=hourly_models,
      misstolerance = 10)
    
    #past_scenarios<-daily_models_past_scenarios
    past_scenarios<-lapply(
      names(daily_models_past_scenarios),
      function(x)
        cbind(daily_models_past_scenarios[[x]],
              hourly_models_past_scenarios[[x]][,names(hourly_models)]))
    names(past_scenarios)<-names(daily_models_past_scenarios)
    
    daily_models_observed<-tempResponse_daily(
      MApfelstetten_temps,
      Start_JDay = Start_JDay,
      End_JDay = End_JDay,
      models=daily_models)
    daily_models_observed<-
      daily_models_observed[which(daily_models_observed$Perc_complete>90),]
    hourly_models_observed<-tempResponse_daily_list(
      MApfelstetten_temps,
      latitude=48.493,
      Start_JDay = Start_JDay,
      End_JDay = End_JDay,
      models=hourly_models,
      misstolerance = 10)
    
    past_observed<-cbind(
      daily_models_observed,
      hourly_models_observed[[1]][,names(hourly_models)])
    
    save_temperature_scenarios(past_scenarios,
                               "data/chill",
                               "MApfelstetten_multichill_historic")
    write.csv(past_observed,
              "data/chill/MApfelstetten_multichill_observed.csv",
              row.names=FALSE)
    
    
    ##same procedure for future data
    RCPs<-c("rcp45","rcp85")
    Times<-c(2050,2090)
    
    for(RCP in RCPs)
      for(Time in Times)
        {
        Temps<-load_temperature_scenarios(
          "data/Weather",
          paste0("MApfelstetten_",Time,"_",RCP))
    
        daily_models_future_scenarios<-tempResponse_list_daily(
          Temps,
          Start_JDay = Start_JDay,
          End_JDay = End_JDay,
          models=daily_models)
        daily_models_future_scenarios<-lapply(
          daily_models_future_scenarios,
          function(x) x[which(x$Perc_complete>90),])
    
        hourly_models_future_scenarios<-
          tempResponse_daily_list(
            Temps,
            latitude=48.493,
            Start_JDay = Start_JDay,
            End_JDay = End_JDay,
            models=hourly_models,
            misstolerance = 10)
    
        future_scenarios<-daily_models_future_scenarios
        future_scenarios<-lapply(
          names(future_scenarios),
          function(x)
            cbind(future_scenarios[[x]],
                  hourly_models_future_scenarios[[x]][,names(hourly_models)]))
        names(future_scenarios)<-names(daily_models_future_scenarios)
    
        chill<-future_scenarios
        save_temperature_scenarios(
          chill,
          "data/chill",
          paste0("MApfelstetten_multichill_",Time,"_",RCP))
    }
    
    
    ##loading and producing chill scenarios
    chill_past_scenarios<-load_temperature_scenarios(
      "data/chill",
      "MApfelstetten_multichill_historic")
    chill_observed<-read_tab("data/chill/MApfelstetten_multichill_observed.csv")
    
    chills <-make_climate_scenario(chill_past_scenarios,
                                   caption = "Historic",
                                   historic_data = chill_observed,
                                   time_series = TRUE)
    
    for(RCP in RCPs)
      for(Time in Times)
        {
        chill<-load_temperature_scenarios(
          "data/chill",
          paste0("MApfelstetten_multichill_",Time,"_",RCP))
        if(RCP=="rcp45") RCPcaption <- "RCP4.5"
        if(RCP=="rcp85") RCPcaption <- "RCP8.5"
        if(Time=="2050") Time_caption <- "2050"
        if(Time=="2090") Time_caption <- "2090"
        chills <-make_climate_scenario(chill,
                                       caption =c(RCPcaption,Time_caption),
                                       add_to = chills)
    }
    
    
    ##heat map for Safe Winter chill
    for(i in 1:length(chills))
       {ch<-chills[[i]]
        if(ch$caption[1]=="Historic")
          {GCMs<-rep("none",length(names(ch$data)))
           RCPs<-rep("none",length(names(ch$data)))
           Years<-as.numeric(ch$labels)
           Scenario<-rep("Historic",length(names(ch$data)))} else
             {GCMs<-names(ch$data)
              RCPs<-rep(ch$caption[1],length(names(ch$data)))
              Years<-rep(as.numeric(ch$caption[2]),length(names(ch$data)))
              Scenario<-rep("Future",length(names(ch$data)))}
        for(nam in names(ch$data))
        {for(met in metrics)
          {temp_res<-data.frame(Metric=met,
                       GCM=GCMs[which(nam==names(ch$data))],
                       RCP=RCPs[which(nam==names(ch$data))],
                       Year=Years[which(nam==names(ch$data))],
                       Result=quantile(ch$data[[nam]][,met],0.1),
                       Scenario=Scenario[which(nam==names(ch$data))])
          if(i==1&nam==names(ch$data)[1]&met==metrics[1])
            results<-temp_res else
              results<-rbind(results,temp_res)
    
        }}}
    
    for(met in metrics)
      results[which(results$Metric==met),"SWC"]<-
        results[which(results$Metric==met),"Result"]/
          results[which(results$Metric==met&results$Year==1980),"Result"]-1
    
    
    
    rng = range(results$SWC)
    
    p_future<-ggplot(results[which(!results$GCM=="none"),],
                     aes(GCM, y=factor(Metric, levels=metrics),
                         fill = SWC)) +
      geom_tile()
    
    p_future
    
    
    p_future <-
      p_future +
      facet_grid(RCP ~ Year) 
    p_future
    
    #make text smaller
    p_future <-
      p_future +
      theme_bw(base_size = 15) +
      theme(axis.text = element_text(size=8))
    
    p_future
    
    library(colorRamps)
    p_future <-
      p_future +
      scale_fill_gradientn(colours=matlab.like(15),
                           labels = scales::percent,
                           limits=rng)
    p_future
    
    
    p_future <-
      p_future  +
      theme(axis.text.x = element_text(angle = 75, hjust = 1, vjust = 1)) +
      labs(fill = "Change in\nSafe Winter Chill\nsince 1975") +
      scale_y_discrete(labels=model_labels) +
      ylab("Chill metric")
    p_future
    
    
    #plot of past scenarios
    
    p_past<-
      ggplot(results[which(results$GCM=="none"),],
             aes(Year, y=factor(Metric, levels=metrics),
                 fill = SWC)) +
      geom_tile()
    
    p_past<-
      p_past +
      theme_bw(base_size = 12) +
      theme(axis.text = element_text(size=8))
    
    p_past<-
      p_past +
      scale_fill_gradientn(colours=matlab.like(15),
                           labels = scales::percent,
                           limits=rng)
    
    p_past<-
      p_past +
      scale_x_continuous(position = "top") 
    
    p_past<-
      p_past +
      labs(fill = "Change in\nSafe Winter Chill\nsince 1975") +
      scale_y_discrete(labels=model_labels) +
      ylab("Chill metric")
    
    
    chill_comp_plot<-
      (p_past +
         p_future +
         plot_layout(guides = "collect",nrow=2, heights=c(1,2))) &
      theme(legend.position = "right",strip.background = element_blank(),
            strip.text = element_text(face = "bold"))
    
    chill_comp_plot

  3. Produce an animated line plot of your results (summarizing Safe Winter Chill across all the GCMs).

    library(chillR)
    ##install.packages("remotes")
    ##install.packages("colorRamps")
    install.packages("gganimate")
    library(remotes)
    remotes::install_github("EduardoFernandezC/dormancyR")
    library(dormancyR)
    library(ggplot2)
    library(reshape2)
    library(colorRamps)
    library(kableExtra)
    library(patchwork)
    library(gganimate)
    
    RCPs<-c("rcp45","rcp85")
    Times<-c(2050,2090)
    
    hourly_models <- list(Chilling_units = chilling_units,
         Low_chill = low_chill_model,
         Modified_Utah = modified_utah_model,
         North_Carolina = north_carolina_model,
         Positive_Utah = positive_utah_model,
         Chilling_Hours = Chilling_Hours,
         Utah_Chill_Units = Utah_Model,
         Chill_Portions = Dynamic_Model)
    daily_models<-list(Rate_of_Chill = rate_of_chill,
                       Chill_Days = chill_days,
                       Exponential_Chill = exponential_chill,
                       Triangula_Chill_Haninnen = triangular_chill_1,
                       Triangular_Chill_Legave = triangular_chill_2)
    
    metrics<-c(names(daily_models),names(hourly_models))
    
    model_labels=c("Rate of Chill",
                   "Chill Days",
                   "Exponential Chill",
                   "Triangular Chill (Häninnen)",
                   "Triangular Chill (Legave)",
                   "Chilling Units",
                   "Low-Chill Chill Units",
                   "Modified Utah Chill Units",
                   "North Carolina Chill Units",
                   "Positive Utah Chill Units",
                   "Chilling Hours",
                   "Utah Chill Units",
                   "Chill Portions")
    
    kable(data.frame(Metric=model_labels,'Function name'=metrics))  %>%
      kable_styling("striped", position = "left",font_size = 10)
    
    MApfelstetten_temps<-read_tab("data/MApfelstetten_temps.csv")
    Temps<-load_temperature_scenarios("data/Weather","MApfelstetten_historic")
    
    
    Start_JDay<-305
    End_JDay<-59
    
    daily_models_past_scenarios<-tempResponse_list_daily(
      Temps,
      Start_JDay = Start_JDay,
      End_JDay = End_JDay,
      models=daily_models)
    
    
    hourly_models_past_scenarios<-tempResponse_daily_list(
      Temps,
      latitude=48.493,
      Start_JDay = Start_JDay,
      End_JDay = End_JDay,
      models=hourly_models,
      misstolerance = 10)
    
    #past_scenarios<-daily_models_past_scenarios
    past_scenarios<-lapply(
      names(daily_models_past_scenarios),
      function(x)
        cbind(daily_models_past_scenarios[[x]],
              hourly_models_past_scenarios[[x]][,names(hourly_models)]))
    names(past_scenarios)<-names(daily_models_past_scenarios)
    
    daily_models_observed<-tempResponse_daily(
      MApfelstetten_temps,
      Start_JDay = Start_JDay,
      End_JDay = End_JDay,
      models=daily_models)
    daily_models_observed<-
      daily_models_observed[which(daily_models_observed$Perc_complete>90),]
    hourly_models_observed<-tempResponse_daily_list(
      MApfelstetten_temps,
      latitude=48.493,
      Start_JDay = Start_JDay,
      End_JDay = End_JDay,
      models=hourly_models,
      misstolerance = 10)
    
    past_observed<-cbind(
      daily_models_observed,
      hourly_models_observed[[1]][,names(hourly_models)])
    
    save_temperature_scenarios(past_scenarios,
                               "data/chill",
                               "MApfelstetten_multichill_historic")
    write.csv(past_observed,
              "data/chill/MApfelstetten_multichill_observed.csv",
              row.names=FALSE)
    
    
    ##same procedure for future data
    RCPs<-c("rcp45","rcp85")
    Times<-c(2050,2090)
    
    for(RCP in RCPs)
      for(Time in Times)
        {
        Temps<-load_temperature_scenarios(
          "data/Weather",
          paste0("MApfelstetten_",Time,"_",RCP))
    
        daily_models_future_scenarios<-tempResponse_list_daily(
          Temps,
          Start_JDay = Start_JDay,
          End_JDay = End_JDay,
          models=daily_models)
        daily_models_future_scenarios<-lapply(
          daily_models_future_scenarios,
          function(x) x[which(x$Perc_complete>90),])
    
        hourly_models_future_scenarios<-
          tempResponse_daily_list(
            Temps,
            latitude=48.493,
            Start_JDay = Start_JDay,
            End_JDay = End_JDay,
            models=hourly_models,
            misstolerance = 10)
    
        future_scenarios<-daily_models_future_scenarios
        future_scenarios<-lapply(
          names(future_scenarios),
          function(x)
            cbind(future_scenarios[[x]],
                  hourly_models_future_scenarios[[x]][,names(hourly_models)]))
        names(future_scenarios)<-names(daily_models_future_scenarios)
    
        chill<-future_scenarios
        save_temperature_scenarios(
          chill,
          "data/chill",
          paste0("MApfelstetten_multichill_",Time,"_",RCP))
    }
    
    
    ##loading and producing chill scenarios
    chill_past_scenarios<-load_temperature_scenarios(
      "data/chill",
      "MApfelstetten_multichill_historic")
    chill_observed<-read_tab("data/chill/MApfelstetten_multichill_observed.csv")
    
    chills <-make_climate_scenario(chill_past_scenarios,
                                   caption = "Historic",
                                   historic_data = chill_observed,
                                   time_series = TRUE)
    
    for(RCP in RCPs)
      for(Time in Times)
        {
        chill<-load_temperature_scenarios(
          "data/chill",
          paste0("MApfelstetten_multichill_",Time,"_",RCP))
        if(RCP=="rcp45") RCPcaption <- "RCP4.5"
        if(RCP=="rcp85") RCPcaption <- "RCP8.5"
        if(Time=="2050") Time_caption <- "2050"
        if(Time=="2090") Time_caption <- "2090"
        chills <-make_climate_scenario(chill,
                                       caption =c(RCPcaption,Time_caption),
                                       add_to = chills)
    }
    
    
    ##heat map for Safe Winter chill
    for(i in 1:length(chills))
       {ch<-chills[[i]]
        if(ch$caption[1]=="Historic")
          {GCMs<-rep("none",length(names(ch$data)))
           RCPs<-rep("none",length(names(ch$data)))
           Years<-as.numeric(ch$labels)
           Scenario<-rep("Historic",length(names(ch$data)))} else
             {GCMs<-names(ch$data)
              RCPs<-rep(ch$caption[1],length(names(ch$data)))
              Years<-rep(as.numeric(ch$caption[2]),length(names(ch$data)))
              Scenario<-rep("Future",length(names(ch$data)))}
        for(nam in names(ch$data))
        {for(met in metrics)
          {temp_res<-data.frame(Metric=met,
                       GCM=GCMs[which(nam==names(ch$data))],
                       RCP=RCPs[which(nam==names(ch$data))],
                       Year=Years[which(nam==names(ch$data))],
                       Result=quantile(ch$data[[nam]][,met],0.1),
                       Scenario=Scenario[which(nam==names(ch$data))])
          if(i==1&nam==names(ch$data)[1]&met==metrics[1])
            results<-temp_res else
              results<-rbind(results,temp_res)
    
        }}}
    
    for(met in metrics)
      results[which(results$Metric==met),"SWC"]<-
        results[which(results$Metric==met),"Result"]/
          results[which(results$Metric==met&results$Year==1980),"Result"]-1
    
    ##animated plot creation
    hist_results<-results[which(results$GCM=="none"),]
    hist_results$RCP<-"RCP4.5"
    hist_results_2<-hist_results
    hist_results_2$RCP<-"RCP8.5"
    hist_results<-rbind(hist_results,hist_results_2)
    
    future_results<-results[which(!results$GCM=="none"),]
    
    GCM_aggregate<-aggregate(
      future_results$SWC,
      by=list(future_results$Metric,future_results$RCP,future_results$Year),
      FUN=mean)
    colnames(GCM_aggregate)<-c("Metric","RCP","Year","SWC")
    
    RCP_Time_series<-rbind(hist_results[,c("Metric","RCP","Year","SWC")],
                           GCM_aggregate)
    
    
    
    chill_change_plot<-
      ggplot(data=RCP_Time_series,
             aes(x=Year,y=SWC,col=factor(Metric,levels=metrics))) +
      geom_line(linewidth=1.3) +
      facet_wrap(~RCP,nrow=2) +
      theme_bw(base_size=12) +
      labs(col = "Change in\nSafe Winter Chill\nsince 1975") +
      scale_color_discrete(labels=model_labels) +
      scale_y_continuous(labels = scales::percent) +
      theme(strip.background = element_blank(),
            strip.text = element_text(face = "bold")) +
      ylab("Safe Winter Chill")
    
    chill_change_plot
    
    
    library(gganimate)
    p<-chill_change_plot + transition_reveal(Year)
    
    ## the following is needed to bypass a writing permission issue - not sure why but
    ## this appears to be working
    
    install.packages("gifski")
    library(gifski)
    animate(p, renderer = gifski_renderer())
    
    anim_save("data/chill_comparison_animation.gif", animation = last_animation())

17 Simple phenology analysis

Exercise:

  1. Provide a brief narrative describing what p-hacking is, and why this is a problematic approach to data analysis.

    In regression fitting, a p-value lower than 0.05 means that the null hypothesis can be rejected (thus there is a significant difference). p-hacking is when researchers handling huge data try out numerous statistical analyses to find and fit patterns or trends that are ‘significant’. These published “cherry-picked” correlations can mislead future research, impede scientific progress, and misguide people by attributing results with the wrong parameters.

  2. Provide a sketch of your causal understanding of the relationship between temperature and bloom dates.

    Bloom dates<-f(chilling, forcing), which are temperature-dependent. Hence, more heating means less chill and more forcing; this leads to earlier bloom. Specifically, due to the increase in temperatures in autumn-winter, the blooming schedules get impaired, a phenomenon supported by the earlier bloom of Japanese pear investigated by Tominaga et al in 2022. The researchers speculate that the lack of chilling (warmer temperatures) caused metabolic and physiological changes leading to earlier bloom dates or ‘flowering disorder’.

  3. What do we need to know to build a process-based model from this?

    From the previous lessons, I learned that process-based models are the most complex models where lots of parameters need to be considered. I also remember that there should be a point at which we draw the line because considering too many factors can multiply errors too; thus we need a moderate number of variables. We need to consider chill and heat accumulation to build the model.

18 Delineating temperature response phases with PLS regression

Exercise:

  1. Briefly explain why you shouldn’t take the results of a PLS regression analysis between temperature and phenology at face value. What do you need in addition in order to make sense of such outputs?

    We must not over-interpret PLS regression analysis (temperature and phenology) results because this cannot prove anything on its own. It is comparable to finding patterns in the walls of the bathroom and interpreting them as meaningful images. To make sense of things, we must interpret trends side by side with concrete, proven, and scientific basis– in this case, the ecological theories. Ecological theories can be used to set expectations, then PLS analysis can be used to support or challenge the ideas.

  2. Replicate the PLS analysis for the Roter Boskoop dataset that you used in a previous lesson.

    library(chillR)
    library(kableExtra)
    Roter<-read_tab("data/Roter_Boskoop_bloom_1958_2019.csv")
    Roter_first<-Roter[,1:2]
    Roter_first[,"Year"]<-substr(Roter_first$First_bloom,1,4)
    Roter_first[,"Month"]<-substr(Roter_first$First_bloom,5,6)
    Roter_first[,"Day"]<-substr(Roter_first$First_bloom,7,8)
    Roter_first<-make_JDay(Roter_first)
    Roter_first<-Roter_first[,c("Pheno_year","JDay")]
    colnames(Roter_first)<-c("Year","pheno")
    
    kable(Roter_first[1:5,]) %>%
      kable_styling("striped", position = "left",font_size = 8)
    
    KA_temps<-read_tab("data/TMaxTMin1958-2019_patched.csv")
    KA_temps<-make_JDay(KA_temps)
    
    kable(KA_temps[1:5,])  %>%
      kable_styling("striped", position = "left",font_size = 8)
    
    
    PLS_results<-PLS_pheno(KA_temps,Roter_first)
    
    kable(PLS_results$PLS_summary[1:5,])  %>%
      kable_styling("striped", position = "left",font_size = 8)
    
    dir.create("pictures")
    plot_PLS(PLS_results,PLS_results_path = "pictures/PLS_out")
    
    
    
    library(ggplot2)
    
    PLS_gg<-PLS_results$PLS_summary
    PLS_gg[,"Month"]<-trunc(PLS_gg$Date/100)
    PLS_gg[,"Day"]<-PLS_gg$Date-PLS_gg$Month*100
    PLS_gg[,"Date"]<-ISOdate(2002,PLS_gg$Month,PLS_gg$Day)
    PLS_gg[which(PLS_gg$JDay<=0),"Date"]<-
      ISOdate(2001,
              PLS_gg$Month[which(PLS_gg$JDay<=0)],
              PLS_gg$Day[which(PLS_gg$JDay<=0)])
    PLS_gg[,"VIP_importance"]<-PLS_gg$VIP>=0.8
    PLS_gg[,"VIP_Coeff"]<-factor(sign(PLS_gg$Coef)*PLS_gg$VIP_importance)
    
    VIP_plot<- ggplot(PLS_gg,aes(x=Date,y=VIP)) +
      geom_bar(stat='identity',aes(fill=VIP_importance))
    VIP_plot <- VIP_plot +
      scale_fill_manual(name="VIP", 
                        labels = c("<0.8", ">0.8"), 
                        values = c("FALSE"="grey", "TRUE"="blue")) +
      theme_bw(base_size=12) +
      theme(axis.text.x = element_blank(),
            axis.ticks.x = element_blank(),
            axis.title.x = element_blank() )
    
    VIP_plot
    
    coeff_plot<- ggplot(PLS_gg,aes(x=Date,y=Coef)) +
      geom_bar(stat='identity',aes(fill=VIP_Coeff)) +
      scale_fill_manual(name="Effect direction", 
                        labels = c("Advancing", "Unimportant","Delaying"), 
                        values = c("-1"="red", "0"="grey","1"="dark green")) +
      theme_bw(base_size=12) +
      ylab("PLS coefficient") +
      theme(axis.text.x = element_blank(),
            axis.ticks.x = element_blank(),
            axis.title.x = element_blank() )
    
    coeff_plot
    
    temp_plot<- ggplot(PLS_gg) +
      geom_ribbon(aes(x=Date,ymin=Tmean-Tstdev,ymax=Tmean+Tstdev),
                  fill="grey") +
      geom_ribbon(aes(x=Date,ymin=Tmean-Tstdev*(VIP_Coeff==-1),
                      ymax=Tmean+Tstdev*(VIP_Coeff==-1)),
                  fill="red") +
      geom_ribbon(aes(x=Date,ymin=Tmean-Tstdev*(VIP_Coeff==1),
                      ymax=Tmean+Tstdev*(VIP_Coeff==1)),
                  fill="dark green") +
      geom_line(aes(x=Date,y=Tmean)) +
      theme_bw(base_size=12) +
      ylab(expression(paste(T[mean]," (°C)")))
    
    temp_plot
    
    library(patchwork)
    
    plot<- (VIP_plot +
              coeff_plot +
              temp_plot +
              plot_layout(ncol=1,
                guides = "collect")
            ) & theme(legend.position = "right",
                      legend.text = element_text(size=8),
                      legend.title = element_text(size=10),
                      axis.title.x=element_blank())
    
    plot
    
    
    ggplot_PLS<-function(PLS_results)
    {
      library(ggplot2)
      PLS_gg<-PLS_results$PLS_summary
      PLS_gg[,"Month"]<-trunc(PLS_gg$Date/100)
      PLS_gg[,"Day"]<-PLS_gg$Date-PLS_gg$Month*100
      PLS_gg[,"Date"]<-ISOdate(2002,PLS_gg$Month,PLS_gg$Day)
      PLS_gg[which(PLS_gg$JDay<=0),"Date"]<-
        ISOdate(2001,
                PLS_gg$Month[which(PLS_gg$JDay<=0)],
                PLS_gg$Day[which(PLS_gg$JDay<=0)])
      PLS_gg[,"VIP_importance"]<-PLS_gg$VIP>=0.8
      PLS_gg[,"VIP_Coeff"]<-factor(sign(PLS_gg$Coef)*PLS_gg$VIP_importance)
    
      VIP_plot<- ggplot(PLS_gg,aes(x=Date,y=VIP)) +
      geom_bar(stat='identity',aes(fill=VIP>0.8)) +
      scale_fill_manual(name="VIP", 
                        labels = c("<0.8", ">0.8"), 
                        values = c("FALSE"="grey", "TRUE"="blue")) +
      theme_bw(base_size=12) +
      theme(axis.text.x = element_blank(),
            axis.ticks.x = element_blank(),
            axis.title.x = element_blank() )
    
      coeff_plot<- ggplot(PLS_gg,aes(x=Date,y=Coef)) +
      geom_bar(stat='identity',aes(fill=VIP_Coeff)) +
      scale_fill_manual(name="Effect direction", 
                        labels = c("Advancing", "Unimportant","Delaying"), 
                        values = c("-1"="red", "0"="grey","1"="dark green")) +
      theme_bw(base_size=12) +
      ylab("PLS coefficient") +
      theme(axis.text.x = element_blank(),
            axis.ticks.x = element_blank(),
            axis.title.x = element_blank() )
    
      temp_plot<- ggplot(PLS_gg) +
        geom_ribbon(aes(x=Date,ymin=Tmean-Tstdev,ymax=Tmean+Tstdev),
                    fill="grey") +
        geom_ribbon(aes(x=Date,ymin=Tmean-Tstdev*(VIP_Coeff==-1),
                        ymax=Tmean+Tstdev*(VIP_Coeff==-1)),
                    fill="red") +
        geom_ribbon(aes(x=Date,ymin=Tmean-Tstdev*(VIP_Coeff==1),
                        ymax=Tmean+Tstdev*(VIP_Coeff==1)),
                    fill="dark green") +
        geom_line(aes(x=Date,y=Tmean)) +
        theme_bw(base_size=12) +
        ylab(expression(paste(T[mean]," (°C)")))
    
      library(patchwork)
      plot<- (VIP_plot +
                coeff_plot +
                temp_plot +
                plot_layout(ncol=1,
                            guides = "collect")
              ) & theme(legend.position = "right",
                        legend.text = element_text(size=8),
                        legend.title = element_text(size=10),
                        axis.title.x=element_blank())
    
    plot
    }
    
    ggplot_PLS(PLS_results)

  3. Write down your thoughts on why we’re not seeing the temperature response pattern we may have expected. What happened to the chill response?

    From the plot, between mid February to April- advancing bloom dates due to high heat accumulation was observed (clear forcing period; almost same trend with the example given for this chapter. Chilling period is not clear, as what was observed in Alexander Lucas pear dataset in Klein-Altendorf. One reason for unclear chilling response is due to cooler climatic conditions–this blurs the relationship between temperature changes and chill accumulation. Klein-Altendorf was one example of a cold location, maybe the Roter-Boskoop data was also captured in a relatively cooler location.

19 Successes and limitations of PLS regression analysis

Exercise:

  1. Briefly explain in what climatic settings we can expect PLS regression to detect the chilling phase - and in what settings this probably won’t work.

    In settings with relatively warm temperatures (in the case of the given examples, Davis), chilling period can be clearly detected by PLS; in relatively cooler settings (e.g. Klein-Alterdorf, Beijing), chill accumulation and temperature relationship is not monotonic, so results can be fuzzy.

  2. How could we overcome this problem?

    Conversion of temperatures to daily rate of chill accumulation. In Chapter 20, this was achieved by using chillR models, particularly with Dynamic Model.

20 PLS regression with agroclimatic metrics

In finding chilling and forcing periods:

-consistent negative coefs

Chilling period-> relationship between bloom dates and daily chill accumulation

Forcing period-> relationship between blooom dates and daily heat accumulation

In delineating, look at results tables, plots, and ecological theory behind it.

Exercise:

  1. Repeat the PLS_chill_force procedure for the ‘Roter Boskoop’ dataset. Include plots of daily chill and heat accumulation.

    library(chillR)
    library(kableExtra)
    
    temps<-read_tab("data/TMaxTMin1958-2019_patched.csv")
    temps_hourly<-stack_hourly_temps(temps,latitude=48.5)
    
    library(dplyr)
    
    temps_hourly<-read_tab("data/TMaxTMin1958-2019_patched.csv") |>
      stack_hourly_temps(latitude=48.5)
    
    head(temps_hourly$hourtemps)
    
    daychill<-daily_chill(hourtemps=temps_hourly,
                          running_mean=1,
                          models = list(Chilling_Hours = Chilling_Hours, Utah_Chill_Units = Utah_Model,
                                        Chill_Portions = Dynamic_Model, GDH = GDH)
    )
    
    head(daychill$daily_chill)
    
    dc<-make_daily_chill_plot2(daychill,metrics=c("Chill_Portions"),cumulative=FALSE,
                               startdate=300,enddate=30,focusyears=c(2008), metriclabels="Chill Portions")
    
    
    dc<-make_daily_chill_plot2(daychill,metrics=c("Chill_Portions"),cumulative=TRUE,
                               startdate=300,enddate=30,focusyears=c(2008), metriclabels="Chill Portions")
    
    
    Roter<-read_tab("data/Roter_Boskoop_bloom_1958_2019.csv")
    Roter_first<-Roter[,1:2]
    Roter_first[,"Year"]<-substr(Roter_first$First_bloom,1,4)
    Roter_first[,"Month"]<-substr(Roter_first$First_bloom,5,6)
    Roter_first[,"Day"]<-substr(Roter_first$First_bloom,7,8)
    Roter_first<-make_JDay(Roter_first)
    Roter_first<-Roter_first[,c("Pheno_year","JDay")]
    colnames(Roter_first)<-c("Year","pheno")
    
    ##extracting PLS results table for CP
    plscf<-PLS_chill_force(daily_chill_obj=daychill,
                           bio_data_frame=Roter_first,
                           split_month=6,
                           chill_models = "Chill_Portions",
                           heat_models = "GDH",
                           runn_means = 11)
    
    plot_PLS(plscf, PLS_results_path= "data/plots/PLS_outputs_11_day")
    
    ##extracting PLS results table for CU
    plscf<-PLS_chill_force(daily_chill_obj=daychill,
                           bio_data_frame=Roter_first,
                           split_month=6,
                           chill_models = "Utah_Chill_Units",
                           heat_models = "GDH",
                           runn_means = 11)
    
    plot_PLS(plscf, PLS_results_path= "data/plots/PLS_outputs_11_day")
    
    ##extracting PLS results table for CH
    plscf<-PLS_chill_force(daily_chill_obj=daychill,
                           bio_data_frame=Roter_first,
                           split_month=6,
                           chill_models = "Chilling_Hours",
                           heat_models = "GDH",
                           runn_means = 11)
    
    plot_PLS(plscf, PLS_results_path= "data/plots/PLS_outputs_11_day")

    Above is the PLS analysis plot (using CP model) applied to the Roter-Boskoop dataset. From here, I can roughly say that chilling period started around early December to February, which is quite similar with the example given in class. Forcing period must be around January to bloom date.

  2. Run PLS_chill_force analyses for all three major chill models. Delineate your best estimates of chilling and forcing phases for all of them.

    library(chillR)
    library(kableExtra)
    library(ggplot2)
    library(dplyr)
    
    temps<-read_tab("data/TMaxTMin1958-2019_patched.csv")
    temps_hourly<-stack_hourly_temps(temps,latitude=48.5)
    
    
    temps_hourly<-read_tab("data/TMaxTMin1958-2019_patched.csv") |>
      stack_hourly_temps(latitude=48.5)
    
    
    Roter<-read_tab("data/Roter_Boskoop_bloom_1958_2019.csv")
    Roter_first<-Roter[,1:2]
    Roter_first[,"Year"]<-substr(Roter_first$First_bloom,1,4)
    Roter_first[,"Month"]<-substr(Roter_first$First_bloom,5,6)
    Roter_first[,"Day"]<-substr(Roter_first$First_bloom,7,8)
    Roter_first<-make_JDay(Roter_first)
    Roter_first<-Roter_first[,c("Pheno_year","JDay")]
    colnames(Roter_first)<-c("Year","pheno")
    
    
    
    daychill<-daily_chill(hourtemps=temps_hourly,
                running_mean=11,
                models = list(Chilling_Hours = Chilling_Hours, Utah_Chill_Units = Utah_Model,
        Chill_Portions = Dynamic_Model, GDH = GDH)
        )
    
    plscf<-PLS_chill_force(daily_chill_obj=daychill,
                           bio_data_frame=Roter_first,
                           split_month=6,
                           chill_models = c("Chilling_Hours", "Utah_Chill_Units", "Chill_Portions"),
                           heat_models = c("GDH"))
    
    ##creation of function for chilling and forcing period shading for CP, CU, and CH
    plot_PLS_chill_force<-function(plscf,
                                   chill_metric="Chill_Portions",
                                   heat_metric="GDH",
                                   chill_label="CP",
                                   heat_label="GDH",
                                   chill_phase=c(-48,62),
                                   heat_phase=c(-5,105.5))
    {
      PLS_gg<-plscf[[chill_metric]][[heat_metric]]$PLS_summary
      PLS_gg[,"Month"]<-trunc(PLS_gg$Date/100)
      PLS_gg[,"Day"]<-PLS_gg$Date-PLS_gg$Month*100
      PLS_gg[,"Date"]<-ISOdate(2002,PLS_gg$Month,PLS_gg$Day)
      PLS_gg[which(PLS_gg$JDay<=0),"Date"]<-ISOdate(2001,PLS_gg$Month[which(PLS_gg$JDay<=0)],PLS_gg$Day[which(PLS_gg$JDay<=0)])
      PLS_gg[,"VIP_importance"]<-PLS_gg$VIP>=0.8
      PLS_gg[,"VIP_Coeff"]<-factor(sign(PLS_gg$Coef)*PLS_gg$VIP_importance)
    
      chill_start_date<-ISOdate(2001,12,31)+chill_phase[1]*24*3600
      chill_end_date<-ISOdate(2001,12,31)+chill_phase[2]*24*3600
      heat_start_date<-ISOdate(2001,12,31)+heat_phase[1]*24*3600
      heat_end_date<-ISOdate(2001,12,31)+heat_phase[2]*24*3600
    
    
      temp_plot<- ggplot(PLS_gg) +
        annotate("rect",
                 xmin = chill_start_date,
                 xmax = chill_end_date,
                 ymin = -Inf,
                 ymax = Inf,
                 alpha = .1,fill = "blue") +
        annotate("rect",
                 xmin = heat_start_date,
                 xmax = heat_end_date,
                 ymin = -Inf,
                 ymax = Inf,
                 alpha = .1,fill = "red") +
        annotate("rect",
                 xmin = ISOdate(2001,12,31) + min(plscf$pheno$pheno,na.rm=TRUE)*24*3600,
                 xmax = ISOdate(2001,12,31) + max(plscf$pheno$pheno,na.rm=TRUE)*24*3600,
                 ymin = -Inf,
                 ymax = Inf,
                 alpha = .1,fill = "black") +
        geom_vline(xintercept = ISOdate(2001,12,31) + median(plscf$pheno$pheno,na.rm=TRUE)*24*3600, linetype = "dashed") +
        geom_ribbon(aes(x=Date,
                        ymin=MetricMean - MetricStdev ,
                        ymax=MetricMean + MetricStdev ),
                    fill="grey") +
        geom_ribbon(aes(x=Date,
                        ymin=MetricMean - MetricStdev * (VIP_Coeff==-1),
                        ymax=MetricMean + MetricStdev * (VIP_Coeff==-1)),
                    fill="red") +
        geom_ribbon(aes(x=Date,
                        ymin=MetricMean - MetricStdev * (VIP_Coeff==1),
                        ymax=MetricMean + MetricStdev * (VIP_Coeff==1)),
                    fill="dark green") +
        geom_line(aes(x=Date,y=MetricMean )) +
        facet_wrap(vars(Type), scales = "free_y",
                   strip.position="left",
                   labeller = labeller(Type = as_labeller(c(Chill=paste0("Chill (",chill_label,")"),Heat=paste0("Heat (",heat_label,")")) )) ) +
        ggtitle("Daily chill and heat accumulation rates") +
        theme_bw(base_size=15) + 
        theme(strip.background = element_blank(),
              strip.placement = "outside",
              strip.text.y = element_text(size =12),
              plot.title = element_text(hjust = 0.5),
              axis.title.y=element_blank()
        )
    
      VIP_plot<- ggplot(PLS_gg,aes(x=Date,y=VIP)) +
        annotate("rect",
                 xmin = chill_start_date,
                 xmax = chill_end_date,
                 ymin = -Inf,
                 ymax = Inf,
                 alpha = .1,fill = "blue") +
        annotate("rect",
                 xmin = heat_start_date,
                 xmax = heat_end_date,
                 ymin = -Inf,
                 ymax = Inf,
                 alpha = .1,fill = "red") +
        annotate("rect",
                 xmin = ISOdate(2001,12,31) + min(plscf$pheno$pheno,na.rm=TRUE)*24*3600,
                 xmax = ISOdate(2001,12,31) + max(plscf$pheno$pheno,na.rm=TRUE)*24*3600,
                 ymin = -Inf,
                 ymax = Inf,
                 alpha = .1,fill = "black") +
        geom_vline(xintercept = ISOdate(2001,12,31) + median(plscf$pheno$pheno,na.rm=TRUE)*24*3600, linetype = "dashed") +
        geom_bar(stat='identity',aes(fill=VIP>0.8)) +
        facet_wrap(vars(Type), scales="free",
                   strip.position="left",
                   labeller = labeller(Type = as_labeller(c(Chill="VIP for chill",Heat="VIP for heat") )) ) +
        scale_y_continuous(limits=c(0,max(plscf[[chill_metric]][[heat_metric]]$PLS_summary$VIP))) +
        ggtitle("Variable Importance in the Projection (VIP) scores") +
        theme_bw(base_size=15) + 
        theme(strip.background = element_blank(),
              strip.placement = "outside",
              strip.text.y = element_text(size =12),
              plot.title = element_text(hjust = 0.5),
              axis.title.y=element_blank()
        ) +
        scale_fill_manual(name="VIP", 
                          labels = c("<0.8", ">0.8"), 
                          values = c("FALSE"="grey", "TRUE"="blue")) +
        theme(axis.text.x = element_blank(),
              axis.ticks.x = element_blank(),
              axis.title.x = element_blank(),
              axis.title.y = element_blank())
    
      coeff_plot<- ggplot(PLS_gg,aes(x=Date,y=Coef)) +
        annotate("rect",
                 xmin = chill_start_date,
                 xmax = chill_end_date,
                 ymin = -Inf,
                 ymax = Inf,
                 alpha = .1,fill = "blue") +
        annotate("rect",
                 xmin = heat_start_date,
                 xmax = heat_end_date,
                 ymin = -Inf,
                 ymax = Inf,
                 alpha = .1,fill = "red") +
        annotate("rect",
                 xmin = ISOdate(2001,12,31) + min(plscf$pheno$pheno,na.rm=TRUE)*24*3600,
                 xmax = ISOdate(2001,12,31) + max(plscf$pheno$pheno,na.rm=TRUE)*24*3600,
                 ymin = -Inf,
                 ymax = Inf,
                 alpha = .1,fill = "black") +
        geom_vline(xintercept = ISOdate(2001,12,31) + median(plscf$pheno$pheno,na.rm=TRUE)*24*3600, linetype = "dashed") +
        geom_bar(stat='identity',aes(fill=VIP_Coeff)) +
        facet_wrap(vars(Type), scales="free",
                   strip.position="left",
                   labeller = labeller(Type = as_labeller(c(Chill="MC for chill",Heat="MC for heat") )) ) +
        scale_y_continuous(limits=c(min(plscf[[chill_metric]][[heat_metric]]$PLS_summary$Coef),
                                    max(plscf[[chill_metric]][[heat_metric]]$PLS_summary$Coef))) +
        ggtitle("Model coefficients (MC)") +
        theme_bw(base_size=15) + 
        theme(strip.background = element_blank(),
              strip.placement = "outside",
              strip.text.y = element_text(size =12),
              plot.title = element_text(hjust = 0.5),
              axis.title.y=element_blank()
        ) +
        scale_fill_manual(name="Effect direction", 
                          labels = c("Advancing", "Unimportant","Delaying"), 
                          values = c("-1"="red", "0"="grey","1"="dark green")) +
        ylab("PLS coefficient") +
        theme(axis.text.x = element_blank(),
              axis.ticks.x = element_blank(),
              axis.title.x = element_blank(),
              axis.title.y = element_blank())
    
      library(patchwork)
    
      plot<- (VIP_plot +
                coeff_plot +
                temp_plot +
                plot_layout(ncol=1,
                            guides = "collect")
      ) & theme(legend.position = "right",
                legend.text = element_text(size=8),
                legend.title = element_text(size=10),
                axis.title.x=element_blank())
    
      plot
    
    }
    
    daychill<-daily_chill(hourtemps=temps_hourly,
                running_mean=11,
                models = list(Chilling_Hours = Chilling_Hours, Utah_Chill_Units = Utah_Model,
        Chill_Portions = Dynamic_Model, GDH = GDH)
        )
    
    plscf<-PLS_chill_force(daily_chill_obj=daychill,
                           bio_data_frame=Roter_first,
                           split_month=6,
                           chill_models = c("Chilling_Hours", "Utah_Chill_Units", "Chill_Portions"),
                           heat_models = c("GDH"))
    
    ##plot (CH)
    plot_PLS_chill_force(plscf,
                         chill_metric="Chilling_Hours",
                         heat_metric="GDH",
                         chill_label="CH",
                         heat_label="GDH",
                         chill_phase=c(0,0),
                         heat_phase=c(0,0))
    ##plot (CU)
    plot_PLS_chill_force(plscf,
                         chill_metric="Utah_Chill_Units",
                         heat_metric="GDH",
                         chill_label="CU",
                         heat_label="GDH",
                         chill_phase=c(0,0),
                         heat_phase=c(0,0))
    
    ##plot (CP)
    plot_PLS_chill_force(plscf,
                         chill_metric="Chill_Portions",
                         heat_metric="GDH",
                         chill_label="CP",
                         heat_label="GDH",
                         chill_phase=c(0,0),
                         heat_phase=c(0,0))
  3. For Chilling Hours model, the chilling period is not clear at all. The model coefficients are unimportant most of the time, so I cannot delineate the chilling period for this. For the forcing period, it is most likely around January to April. I consulted the table of results for the dates, and I arrived at JDay 8 (January 8) to JDay 105.5 (April 15).PLS Chillforce (Chilling Hours)

    For Utah Chilling Units, chilling period is most likely around JDays -15 to 66. Forcing oeriod is most likely from JDays 10 to 105.5.

    For Chilling Portions model, chilling period is most likely around JDays -25 to 62, although I must say that the relationship is not that strong here. Forcing period is most likely from JDays 8 to 105.5; this, in contrast with the former, is very easy to spot.

  4. Plot results for all three analyses, including shaded plot areas for the chilling and forcing periods you estimated.

    library(chillR)
    library(kableExtra)
    library(ggplot2)
    library(dplyr)
    
    temps<-read_tab("data/TMaxTMin1958-2019_patched.csv")
    temps_hourly<-stack_hourly_temps(temps,latitude=48.5)
    
    
    temps_hourly<-read_tab("data/TMaxTMin1958-2019_patched.csv") |>
      stack_hourly_temps(latitude=48.5)
    
    
    Roter<-read_tab("data/Roter_Boskoop_bloom_1958_2019.csv")
    Roter_first<-Roter[,1:2]
    Roter_first[,"Year"]<-substr(Roter_first$First_bloom,1,4)
    Roter_first[,"Month"]<-substr(Roter_first$First_bloom,5,6)
    Roter_first[,"Day"]<-substr(Roter_first$First_bloom,7,8)
    Roter_first<-make_JDay(Roter_first)
    Roter_first<-Roter_first[,c("Pheno_year","JDay")]
    colnames(Roter_first)<-c("Year","pheno")
    
    
    
    daychill<-daily_chill(hourtemps=temps_hourly,
                running_mean=11,
                models = list(Chilling_Hours = Chilling_Hours, Utah_Chill_Units = Utah_Model,
        Chill_Portions = Dynamic_Model, GDH = GDH)
        )
    
    plscf<-PLS_chill_force(daily_chill_obj=daychill,
                           bio_data_frame=Roter_first,
                           split_month=6,
                           chill_models = c("Chilling_Hours", "Utah_Chill_Units", "Chill_Portions"),
                           heat_models = c("GDH"))
    
    ##creation of function for chilling and forcing period shading for CP, CU, and CH
    plot_PLS_chill_force<-function(plscf,
                                   chill_metric="Chill_Portions",
                                   heat_metric="GDH",
                                   chill_label="CP",
                                   heat_label="GDH",
                                   chill_phase=c(-48,62),
                                   heat_phase=c(-5,105.5))
    {
      PLS_gg<-plscf[[chill_metric]][[heat_metric]]$PLS_summary
      PLS_gg[,"Month"]<-trunc(PLS_gg$Date/100)
      PLS_gg[,"Day"]<-PLS_gg$Date-PLS_gg$Month*100
      PLS_gg[,"Date"]<-ISOdate(2002,PLS_gg$Month,PLS_gg$Day)
      PLS_gg[which(PLS_gg$JDay<=0),"Date"]<-ISOdate(2001,PLS_gg$Month[which(PLS_gg$JDay<=0)],PLS_gg$Day[which(PLS_gg$JDay<=0)])
      PLS_gg[,"VIP_importance"]<-PLS_gg$VIP>=0.8
      PLS_gg[,"VIP_Coeff"]<-factor(sign(PLS_gg$Coef)*PLS_gg$VIP_importance)
    
      chill_start_date<-ISOdate(2001,12,31)+chill_phase[1]*24*3600
      chill_end_date<-ISOdate(2001,12,31)+chill_phase[2]*24*3600
      heat_start_date<-ISOdate(2001,12,31)+heat_phase[1]*24*3600
      heat_end_date<-ISOdate(2001,12,31)+heat_phase[2]*24*3600
    
    
      temp_plot<- ggplot(PLS_gg) +
        annotate("rect",
                 xmin = chill_start_date,
                 xmax = chill_end_date,
                 ymin = -Inf,
                 ymax = Inf,
                 alpha = .1,fill = "blue") +
        annotate("rect",
                 xmin = heat_start_date,
                 xmax = heat_end_date,
                 ymin = -Inf,
                 ymax = Inf,
                 alpha = .1,fill = "red") +
        annotate("rect",
                 xmin = ISOdate(2001,12,31) + min(plscf$pheno$pheno,na.rm=TRUE)*24*3600,
                 xmax = ISOdate(2001,12,31) + max(plscf$pheno$pheno,na.rm=TRUE)*24*3600,
                 ymin = -Inf,
                 ymax = Inf,
                 alpha = .1,fill = "black") +
        geom_vline(xintercept = ISOdate(2001,12,31) + median(plscf$pheno$pheno,na.rm=TRUE)*24*3600, linetype = "dashed") +
        geom_ribbon(aes(x=Date,
                        ymin=MetricMean - MetricStdev ,
                        ymax=MetricMean + MetricStdev ),
                    fill="grey") +
        geom_ribbon(aes(x=Date,
                        ymin=MetricMean - MetricStdev * (VIP_Coeff==-1),
                        ymax=MetricMean + MetricStdev * (VIP_Coeff==-1)),
                    fill="red") +
        geom_ribbon(aes(x=Date,
                        ymin=MetricMean - MetricStdev * (VIP_Coeff==1),
                        ymax=MetricMean + MetricStdev * (VIP_Coeff==1)),
                    fill="dark green") +
        geom_line(aes(x=Date,y=MetricMean )) +
        facet_wrap(vars(Type), scales = "free_y",
                   strip.position="left",
                   labeller = labeller(Type = as_labeller(c(Chill=paste0("Chill (",chill_label,")"),Heat=paste0("Heat (",heat_label,")")) )) ) +
        ggtitle("Daily chill and heat accumulation rates") +
        theme_bw(base_size=15) + 
        theme(strip.background = element_blank(),
              strip.placement = "outside",
              strip.text.y = element_text(size =12),
              plot.title = element_text(hjust = 0.5),
              axis.title.y=element_blank()
        )
    
      VIP_plot<- ggplot(PLS_gg,aes(x=Date,y=VIP)) +
        annotate("rect",
                 xmin = chill_start_date,
                 xmax = chill_end_date,
                 ymin = -Inf,
                 ymax = Inf,
                 alpha = .1,fill = "blue") +
        annotate("rect",
                 xmin = heat_start_date,
                 xmax = heat_end_date,
                 ymin = -Inf,
                 ymax = Inf,
                 alpha = .1,fill = "red") +
        annotate("rect",
                 xmin = ISOdate(2001,12,31) + min(plscf$pheno$pheno,na.rm=TRUE)*24*3600,
                 xmax = ISOdate(2001,12,31) + max(plscf$pheno$pheno,na.rm=TRUE)*24*3600,
                 ymin = -Inf,
                 ymax = Inf,
                 alpha = .1,fill = "black") +
        geom_vline(xintercept = ISOdate(2001,12,31) + median(plscf$pheno$pheno,na.rm=TRUE)*24*3600, linetype = "dashed") +
        geom_bar(stat='identity',aes(fill=VIP>0.8)) +
        facet_wrap(vars(Type), scales="free",
                   strip.position="left",
                   labeller = labeller(Type = as_labeller(c(Chill="VIP for chill",Heat="VIP for heat") )) ) +
        scale_y_continuous(limits=c(0,max(plscf[[chill_metric]][[heat_metric]]$PLS_summary$VIP))) +
        ggtitle("Variable Importance in the Projection (VIP) scores") +
        theme_bw(base_size=15) + 
        theme(strip.background = element_blank(),
              strip.placement = "outside",
              strip.text.y = element_text(size =12),
              plot.title = element_text(hjust = 0.5),
              axis.title.y=element_blank()
        ) +
        scale_fill_manual(name="VIP", 
                          labels = c("<0.8", ">0.8"), 
                          values = c("FALSE"="grey", "TRUE"="blue")) +
        theme(axis.text.x = element_blank(),
              axis.ticks.x = element_blank(),
              axis.title.x = element_blank(),
              axis.title.y = element_blank())
    
      coeff_plot<- ggplot(PLS_gg,aes(x=Date,y=Coef)) +
        annotate("rect",
                 xmin = chill_start_date,
                 xmax = chill_end_date,
                 ymin = -Inf,
                 ymax = Inf,
                 alpha = .1,fill = "blue") +
        annotate("rect",
                 xmin = heat_start_date,
                 xmax = heat_end_date,
                 ymin = -Inf,
                 ymax = Inf,
                 alpha = .1,fill = "red") +
        annotate("rect",
                 xmin = ISOdate(2001,12,31) + min(plscf$pheno$pheno,na.rm=TRUE)*24*3600,
                 xmax = ISOdate(2001,12,31) + max(plscf$pheno$pheno,na.rm=TRUE)*24*3600,
                 ymin = -Inf,
                 ymax = Inf,
                 alpha = .1,fill = "black") +
        geom_vline(xintercept = ISOdate(2001,12,31) + median(plscf$pheno$pheno,na.rm=TRUE)*24*3600, linetype = "dashed") +
        geom_bar(stat='identity',aes(fill=VIP_Coeff)) +
        facet_wrap(vars(Type), scales="free",
                   strip.position="left",
                   labeller = labeller(Type = as_labeller(c(Chill="MC for chill",Heat="MC for heat") )) ) +
        scale_y_continuous(limits=c(min(plscf[[chill_metric]][[heat_metric]]$PLS_summary$Coef),
                                    max(plscf[[chill_metric]][[heat_metric]]$PLS_summary$Coef))) +
        ggtitle("Model coefficients (MC)") +
        theme_bw(base_size=15) + 
        theme(strip.background = element_blank(),
              strip.placement = "outside",
              strip.text.y = element_text(size =12),
              plot.title = element_text(hjust = 0.5),
              axis.title.y=element_blank()
        ) +
        scale_fill_manual(name="Effect direction", 
                          labels = c("Advancing", "Unimportant","Delaying"), 
                          values = c("-1"="red", "0"="grey","1"="dark green")) +
        ylab("PLS coefficient") +
        theme(axis.text.x = element_blank(),
              axis.ticks.x = element_blank(),
              axis.title.x = element_blank(),
              axis.title.y = element_blank())
    
      library(patchwork)
    
      plot<- (VIP_plot +
                coeff_plot +
                temp_plot +
                plot_layout(ncol=1,
                            guides = "collect")
      ) & theme(legend.position = "right",
                legend.text = element_text(size=8),
                legend.title = element_text(size=10),
                axis.title.x=element_blank())
    
      plot
    
    }
    
    daychill<-daily_chill(hourtemps=temps_hourly,
                running_mean=11,
                models = list(Chilling_Hours = Chilling_Hours, Utah_Chill_Units = Utah_Model,
        Chill_Portions = Dynamic_Model, GDH = GDH)
        )
    
    plscf<-PLS_chill_force(daily_chill_obj=daychill,
                           bio_data_frame=Roter_first,
                           split_month=6,
                           chill_models = c("Chilling_Hours", "Utah_Chill_Units", "Chill_Portions"),
                           heat_models = c("GDH"))
    ##shaded plot (CH)
    plot_PLS_chill_force(plscf,
                         chill_metric="Chilling_Hours",
                         heat_metric="GDH",
                         chill_label="CH",
                         heat_label="GDH",
                         chill_phase=c(0,0),
                         heat_phase=c(8,105.5))
    
    
    ##shaded plot (CU)
    plot_PLS_chill_force(plscf,
                         chill_metric="Utah_Chill_Units",
                         heat_metric="GDH",
                         chill_label="CU",
                         heat_label="GDH",
                         chill_phase=c(-15,66),
                         heat_phase=c(10,105.5))
    
    ##shaded plot (CP)
    plot_PLS_chill_force(plscf,
                         chill_metric="Chill_Portions",
                         heat_metric="GDH",
                         chill_label="CP",
                         heat_label="GDH",
                         chill_phase=c(-25,62),
                         heat_phase=c(8,105.5))

    For CH model:

    For CU model:

    For CP model:

21 Examples of PLS regression with agroclimatic metrics

Exercise:

  1. Look across all the PLS results presented above. Can you detect a pattern in where chilling and forcing periods could be delineated clearly, and where this attempt failed?

    In cold conditions (Beijing), it is very hard to delineate chilling period from the graphs, but forcing period is clear.

    In slightly warmer conditions like Klein-Altendorf, it was still difficult to see the chilling period still.

    In warmer conditions, delineating chilling period may be easy due to very strong relationships (high VIP and consistently very low model coefficients), they seem to not agree with expected range of dates (extends beyond normal). As conditions get warmer, chilling periods get easier to spot and get longer, while forcing gets harder to delineate.

    Based on this, PLS may not be that useful in extreme conditions (very cold and very warm). We saw that the graphs gave confusing trends and sometimes outrageous results compared to expected. Chilling periods are hard to see in colder conditions, and forcing periods are hard to delineate in warmer settings.

  2. Can you think about possible reasons for the success or failure of PLS analysis based on agroclimatic metrics (if so, write down your thoughts)?

    In moderate climates (adequate heat and chill during the right time), the chilling and forcing pperiods are easy to delineate using PLS analysis. This might be because they just get enough chill and hit to fulfill in their development stages, so the trends look clear (hits dates we expect based on theory). I think extreme cold and warm settings tend to make things uninterpretable for PLS analysis because there are many possible significant processes, effects ,and relationships arising which are not considered in the models.

22 Why PLS doesn’t always work

Exercise:

  1. Produce chill and heat model sensitivity plots for the location you focused on in previous exercises.

    library(chillR)
    library(dplyr)
    library(lubridate)
    library(ggplot2)
    library(colorRamps)
    library(lubridate) # to get convenient date functionality
    library(purrr) # to get the pluck function
    
    
    ##function for chill model sensitivity
    Chill_model_sensitivity<-
      function(latitude,
               temp_models=list(Dynamic_Model=Dynamic_Model,GDH=GDH),
               month_range=c(10,11,12,1,2,3),
               Tmins=c(-10:20),
               Tmaxs=c(-5:30))
      {
      metrics<-as.list(rep(NA,length(temp_models)))
      names(metrics)<-names(temp_models)
    
      first_run<-TRUE
      for(mon in month_range)
      {weather<-make_all_day_table(
        data.frame(Year=c(2001,2001),
                   Month=c(mon,mon),
                   Day=c(1,days_in_month(mon)),
                   Tmin=c(0,0),
                   Tmax=c(0,0)))
      for(tmin in Tmins)
        for(tmax in Tmaxs)
          if(tmax>=tmin)
          {
            hourtemps <- weather %>%
              mutate(Tmin=tmin, Tmax=tmax) %>% 
              stack_hourly_temps(latitude=latitude) %>% 
              pluck(1,"Temp")         
    
               metric <- lapply(temp_models,
                               function(x) do.call(x,list(hourtemps)) %>% 
                                 tail(1)/(length(hourtemps)/24))
    
              if(first_run)
                {for(tm in 1:length(temp_models))
                  metrics[[tm]]<-metric[[tm]]
                mins<-tmin
                maxs<-tmax
                month<-mon
                first_run<-FALSE} else
                  {for(tm in 1:length(temp_models))
                    metrics[[tm]]<-append(metrics[[tm]],metric[[tm]])
                  mins<-append(mins,tmin)
                  maxs<-append(maxs,tmax)
                  month<-append(month,mon)}
          }
      }
      results<-cbind(data.frame(Month=month,Tmin=mins,Tmax=maxs),
                     as.data.frame(metrics))
      results
    }
    
    
    Chill_sensitivity_temps<-function(chill_model_sensitivity_table,
                                      temperatures,
                                      temp_model,
                                      month_range=c(10,11,12,1,2,3),
                                      Tmins=c(-10:20),
                                      Tmaxs=c(-5:30),
                                      legend_label="Chill/day (CP)")
    {
      library(ggplot2)
      library(colorRamps)
    
      cmst<-filter(chill_model_sensitivity_table,Month %in% month_range) %>% 
        mutate(Month_names = factor(Month,
                                    levels=month_range,
                                    labels=month.name[month_range]))
    
    
      DM_sensitivity<-
        ggplot(cmst,
               aes(x=Tmin,y=Tmax,fill=.data[[temp_model]])) +
        geom_tile() +
        scale_fill_gradientn(colours=alpha(matlab.like(15), alpha = .5),
                             name=legend_label) +
        coord_cartesian(xlim=range(Tmins),
                        ylim=range(Tmaxs))
    
      temperatures<-temperatures %>% 
        filter(Month %in% month_range) %>% 
        filter(Tmax>=Tmin) %>% 
        mutate(Month_names=factor(Month,
                                  levels=month_range,
                                  labels=month.name[month_range])  )
    
    
      DM_sensitivity +
        geom_point(data=temperatures,
                   aes(x=Tmin,y=Tmax,fill=NULL,color="Temperature"),
                   size=0.2) +
        facet_wrap(vars(Month_names)) +
        scale_color_manual(values = "black",
                           labels = "Daily temperature \nextremes (°C)",
                           name="Observed at site" ) +
        guides(fill = guide_colorbar(order = 1),
               color = guide_legend(order = 2)) +
        ylab("Tmax (°C)") +
        xlab("Tmin (°C)") + 
        theme_bw(base_size=12)
    
    }
    
    
    
    Model_sensitivities_MApfelstetten<-
      Chill_model_sensitivity(latitude=48.5,
                              temp_models=list(Dynamic_Model=Dynamic_Model,
                                               GDH=GDH),
                              month_range=c(10:12,1:5))
    
    write.csv(Model_sensitivities_MApfelstetten,
              "data/Model_sensitivities_MApfelstetten.csv",row.names = FALSE)
    
    
    Model_sensitivities_MApfelstetten<-read.csv("data/Model_sensitivities_MApfelstetten.csv")
    
    
    
    MApfelstetten_weather<-read_tab("data/MApfelstetten_weather.csv")
    
    
    
    Chill_sensitivity_temps(Model_sensitivities_MApfelstetten,
                            MApfelstetten_weather,
                            temp_model="Dynamic_Model",
                            month_range=c(10,11,12,1,2,3),
                            legend_label="Chill per day \n(Chill Portions)") +
      ggtitle("Chill model sensitivity near Münsingen Apfelstetten, Germany")
    
    
    
    
    Chill_sensitivity_temps(Model_sensitivities_MApfelstetten,
                            MApfelstetten_weather,
                            temp_model="GDH",
                            month_range=c(12,1:5),
                            legend_label="Heat per day \n(GDH)") +
      ggtitle("Heat model sensitivity near Münsingen Apfelstetten, Germany")

    January and February show the most variation which makes it suitable for PLS analysis, but all the other months seem to have most days within the optimal chilling conditions.

    The GDH model sensitivity plot above shows good variation which makes it clear why it was very easy to spot the forcing period in the PLS analysis.

  2. Look at the plots for the agroclimate-based PLS analyses of the ‘Alexander Lucas’ and ‘Roter Boskoop’ datasets. Provide your best estimates of the chilling and forcing phases.

    Pasting here the PLS Analysis plot for ‘Alexander Lucas’ using the CP model (from the Tree Phenology using R module Chapter 20.4):

    I also looked at the PLS analysis results here, and I think the chilling period is around JDay -31 to JDay 59. Forcing period is most probably from JDay 1 to 105.5 (bloom date).

    Pasting here the shaded PLS Analysis plot for ‘Roter Boskoop’ using the CP model from my work for Chapter 20:

    Chilling period is most likely around JDays -25 to 62, although I must say that the relationship is not that strong here. Forcing period is most likely from JDays 8 to 105.5; this, in contrast with the former, is very easy to spot.

23 Evaluating PLS outputs

Exercise:

  1. Reproduce the analysis for the ‘Roter Boskoop’ dataset.

    library(chillR)
    library(dplyr)
    library(purrr)
    
    
    Roter_first<-read_tab("data/Roter_Boskoop_bloom_1958_2019.csv")[,1:2] %>% 
      mutate(Year=substr(First_bloom,1,4),
             Month=substr(First_bloom,5,6),
             Day=substr(First_bloom,7,8)) %>% 
      make_JDay() %>%
      select(Pheno_year,JDay) %>%
      rename(Year=1,pheno=2)
    
    
    temps<-read_tab("data/TMaxTMin1958-2019_patched.csv")
    
    temps_hourly<-read_tab("data/TMaxTMin1958-2019_patched.csv") %>%
      stack_hourly_temps(latitude=50.6)
    
    daychill<-daily_chill(hourtemps=temps_hourly,
                running_mean=1,
                models = list(Chilling_Hours = Chilling_Hours,
                              Utah_Chill_Units = Utah_Model,
                              Chill_Portions = Dynamic_Model,
                              GDH = GDH)
        )
    
    
    library(ggplot2)
    
    
    chill_phase<-c(340,62)
    heat_phase<-c(357,105.5)
    
    chill<-tempResponse(hourtemps = temps_hourly,
                        Start_JDay = chill_phase[1],
                        End_JDay = chill_phase[2],
                        models = list(Chill_Portions = Dynamic_Model),
                        misstolerance=10)
    
    heat<-tempResponse(hourtemps = temps_hourly,
                        Start_JDay = heat_phase[1],
                        End_JDay = heat_phase[2],
                        models = list(GDH = GDH))
    
    ##plotting 
    ggplot(data=chill,aes(x=Chill_Portions)) +
      geom_histogram() +
      ggtitle("Chill accumulation during endodormancy (Chill Portions)") +
      xlab("Chill accumulation (Chill Portions)") +
      ylab("Frequency between 1958 and 2019") +
      theme_bw(base_size = 12)
    
    ggplot(data=heat,aes(x=GDH)) +
      geom_histogram() +
      ggtitle("Heat accumulation during ecodormancy (GDH)") +
      xlab("Heat accumulation (Growing Degree Hours)") +
      ylab("Frequency between 1958 and 2019") +
      theme_bw(base_size = 12)
    
    
    chill_requirement <- mean(chill$Chill_Portions) ##[1] 56.16143
    chill_req_error <- sd(chill$Chill_Portions)  ##[1] 7.126878
    
    
    median(chill$Chill_Portions) ##[1] 57.6239
    
    heat_requirement <- mean(heat$GDH)   ##[1] 3564.06
    
    
    heat_req_error <- sd(heat$GDH)  ##[1] 1428.505
    
    
    ##function creation for the phenological trend
    pheno_trend_ggplot<-function(temps,
                                 pheno,
                                 chill_phase,
                                 heat_phase,
                                 phenology_stage="Bloom")
    {
      library(fields)
      library(reshape2)
      library(metR)
      library(ggplot2)
      library(colorRamps)
    
      # first, a sub-function (function defined within a function) to
      # compute the temperature means
    
      mean_temp_period<-function(temps,
                                 start_JDay,
                                 end_JDay, 
                                 end_season = end_JDay)
        { temps_JDay<-make_JDay(temps)
        temps_JDay[,"Season"]<-temps_JDay$Year
        if(start_JDay>end_season)
          temps_JDay$Season[which(temps_JDay$JDay>=start_JDay)]<-
          temps_JDay$Year[which(temps_JDay$JDay>=start_JDay)]+1
        if(start_JDay>end_season)
          sub_temps<-subset(temps_JDay,JDay<=end_JDay|JDay>=start_JDay)
        if(start_JDay<=end_JDay)
          sub_temps<-subset(temps_JDay,JDay<=end_JDay&JDay>=start_JDay)
        mean_temps<-aggregate(sub_temps[,c("Tmin","Tmax")],
                              by=list(sub_temps$Season),
                              FUN=function(x) mean(x, na.rm=TRUE))
        mean_temps[,"n_days"]<-aggregate(sub_temps[,"Tmin"],
                                         by=list(sub_temps$Season),
                                         FUN=length)[,2]
        mean_temps[,"Tmean"]<-(mean_temps$Tmin+mean_temps$Tmax)/2
        mean_temps<-mean_temps[,c(1,4,2,3,5)]
        colnames(mean_temps)[1]<-"End_year"
        return(mean_temps)
        }
    
      mean_temp_chill<-mean_temp_period(temps = temps,
                                        start_JDay = chill_phase[1],
                                        end_JDay = chill_phase[2],
                                        end_season = heat_phase[2])
    
      mean_temp_heat<-mean_temp_period(temps = temps,
                                       start_JDay = heat_phase[1],
                                       end_JDay = heat_phase[2],
                                       end_season = heat_phase[2])
    
      mean_temp_chill<-
        mean_temp_chill[which(mean_temp_chill$n_days >= 
                                max(mean_temp_chill$n_days)-1),]
      mean_temp_heat<-
        mean_temp_heat[which(mean_temp_heat$n_days >= 
                               max(mean_temp_heat$n_days)-1),]
      mean_chill<-mean_temp_chill[,c("End_year","Tmean")]
      colnames(mean_chill)[2]<-"Tmean_chill"
      mean_heat<-mean_temp_heat[,c("End_year","Tmean")]
      colnames(mean_heat)[2]<-"Tmean_heat"
      phase_Tmeans<-merge(mean_chill,mean_heat, by="End_year")
    
      colnames(pheno)<-c("End_year","pheno")
      Tmeans_pheno<-merge(phase_Tmeans,pheno, by="End_year")
    
      # Kriging interpolation
      k<-Krig(x=as.matrix(Tmeans_pheno[,c("Tmean_chill","Tmean_heat")]),
              Y=Tmeans_pheno$pheno)
      pred<-predictSurface(k)
      colnames(pred$z)<-pred$y
      rownames(pred$z)<-pred$x
      melted<-melt(pred$z)
      colnames(melted)<-c("Tmean_chill","Tmean_heat","value")
    
      ggplot(melted,aes(x=Tmean_chill,y=Tmean_heat,z=value)) +
        geom_contour_fill(bins=60) +
        scale_fill_gradientn(colours=alpha(matlab.like(15)),
                             name=paste(phenology_stage,"date \n(day of the year)")) +
        geom_contour(col="black") +
        geom_text_contour(stroke = 0.2) +
        geom_point(data=Tmeans_pheno,
                   aes(x=Tmean_chill,y=Tmean_heat,z=NULL),
                   size=0.7)  +
        ylab(expression("Forcing phase "*T[mean]*" ("*degree*C*")")) +
        xlab(expression("Chilling phase "*T[mean]*" ("*degree*C*")")) +
        theme_bw(base_size=15)
    }
    
    
    chill_phase<-c(340,62)
    heat_phase<-c(357,105.5)
    
    pheno_trend_ggplot(temps=temps,
                       pheno=Roter_first,
                       chill_phase=chill_phase,
                       heat_phase=heat_phase,
                       phenology_stage="Bloom")

    I will not repeat the creation of the shaded PLS analysis plot I did for ‘Roter Boskoop’ dataset, I just paste it again here:

    From the analysis, the following values were obtained: chilling requirement (chill_requirement) is around 56.2 CP with an estimated error (chill_req_error) of 7.1 CP. The heat requirement (heat_requirement) is approximately 3564 GDH, with an estimated error (heat_req_error) of 1428 GDH. In temperate regions, chilling requirements range from 64.7 to 79.4 CP, and heat requirements are from 3292.0–8383.6 GDH (Xu et al., 2020); which means the Roter Boskoop’s chilling requirement is lower than normally required by temperate trees, while its heat requirement is just a bit higher than the lower bounds. The phenological plot showed very narrow diagonal shape (strong temperature correlations during eco and endo-dormancy stages), this one is more narrow compared to the ‘Alexander Lucas’ plot. Warm forcing phase led to speeding up of bloom (as early as JDay 100), and cool forcing period led to bloom delay of up to 20 days.

  2. We’ve looked at data from a number of locations so far. How would you expect this surface plot to look like in Beijing? And how should it look in Tunisia?

    Beijing has very cold winters, so probably it would just be an extreme, more pronounced version of the thin, diagonal Klein-Altendorf plot. I would be expecting that cool forcing phase would delay bloom (forcing effect is more dominant in cold conditions), and with warm forcing phase comes earlier bloom. Blue would show up in the upper right side and reds will be on the lower left side.

    Tunisia is quite warm, and one example already given with relatively warm conditions (e.g. Davis), and it had very broad surface plot. Maybe Tunisia’s would be a bit more extremes since it is hotter. A warm chilling phase with cool forcing would delay bloom dates (red on lower right side). The opposite would lead to earlier bloom (blue at top left).

24 The relative importance of chill and heat

Exercise:

Describe the temperature response hypothesis outlined in this chapter.

In ecodormancy stage, bud breaks when target heat accumulation is reached. Endodormancy release happens when target chill requirement is reached. The temperature response theory suggests that shorter dormancy arises from exposure to enough chill during endodormancy stage and enough heat during ecodormancy. Summary:

26 Making valid tree phenology models

Exercise:

  1. Explain the difference between output validation and process validation.

    Output validation only compares (RMSE, Model Efficiency, etc.) the outputs of the models vs observed values; hence neglecting the process (the how’s and why’s). Process validation meshes models with the relevant processes involved; in the context of tree phenology models, the phenological processes should be the reference when doing process validation.

  2. Explain what a validity domain is and why it is important to consider this whenever we want to use our model to forecast something.

    Validity domain is the specific set or range of conditions the model can be used for with acceptable degree of correctness in its predictions.

  3. What is validation for purpose?

    It is the evaluation of the aim and use cases of the model prior to creating it.

  4. How can we ensure that our model is suitable for the predictions we want to make?

    By checking first if the conditions are within the validity domain of the model. If the conditions are not met, then we should not apply the model as is.

27 The PhenoFlex model

Exercise:

  1. Parameterize the PhenoFlex model for `Roter Boskoop’ apples.

    I encountered an error here with the ‘phenologyFitter’ function. Apparently it checks that the seasons list is the same number with the number of years in the bloom dates dataset. And upon inspection, ‘Roter Boskoop’ only had data til 2018, so I had to change that in the ‘SeasonsList’.

    library(chillR)
    library(ggplot2)
    library(kableExtra)
    
    CKA_weather<-read_tab("data/TMaxTMin1958-2019_patched.csv")
    Roter<-read_tab("data/Roter_Boskoop_bloom_1958_2019.csv")
    Roter_first<-Roter[,1:2]
    Roter_first[,"Year"]<-substr(Roter_first$First_bloom,1,4)
    Roter_first[,"Month"]<-substr(Roter_first$First_bloom,5,6)
    Roter_first[,"Day"]<-substr(Roter_first$First_bloom,7,8)
    Roter_first<-make_JDay(Roter_first)
    Roter_first<-Roter_first[,c("Pheno_year","JDay")]
    colnames(Roter_first)<-c("Year","pheno")
    hourtemps <- stack_hourly_temps(CKA_weather, latitude=50.6)
    
    
    # here's the order of the parameters (from the helpfile of the
    # PhenoFlex_GDHwrapper function)
    #          yc,  zc,  s1, Tu,    E0,      E1,     A0,         A1,   Tf, Tc, Tb,  slope
    par <-   c(40, 190, 0.5, 25, 3372.8,  9900.3, 6319.5, 5.939917e13,  4, 36,  4,  1.60)
    upper <- c(41, 200, 1.0, 30, 4000.0, 10000.0, 7000.0,       6.e13, 10, 40, 10, 50.00)
    lower <- c(38, 180, 0.1, 0 , 3000.0,  9000.0, 6000.0,       5.e13,  0,  0,  0,  0.05)
    
    
    
    SeasonList <- genSeasonList(hourtemps$hourtemps, mrange = c(8, 6), years=c(1959:2018))
    
    Fit_res <- phenologyFitter(par.guess=par,
                               modelfn = PhenoFlex_GDHwrapper,
                               bloomJDays=Roter_first$pheno[which(Roter_first$Year>1958)],
                               SeasonList=SeasonList,
                               lower=lower,
                               upper=upper,
                               control=list(smooth=FALSE, verbose=FALSE, maxit=1000,
                                            nb.stop.improvement=5))
    
    
    Roter_par<-Fit_res$par
    
    write.csv(Roter_par,"data/PhenoFlex_parameters_Roter_Boskoop.csv")

    V1

    V2

    NA

    ,“x”

    1

    ,40.5138147244476

    2

    ,199.722028620816

    3

    ,0.18019167088005

    4

    ,24.9705987499286

    5

    ,3372.29197129426

  2. Produce plots of predicted vs. observed bloom dates and distribution of prediction errors.

    library(ggplot2)
    library(tidyr)
    require(reshape2)
    library(chillR)
    
    Roter_par<-read_tab("data/PhenoFlex_parameters_Roter_Boskoop.csv")[,2]
    
    SeasonList <- genSeasonList(hourtemps$hourtemps, mrange = c(8, 6), years=c(1959:2018))
    
    Roter_PhenoFlex_predictions<-Roter_first[which(Roter_first$Year>1958),]
    
    for(y in 1:length(Roter_PhenoFlex_predictions$Year))
       Roter_PhenoFlex_predictions$predicted[y]<-PhenoFlex_GDHwrapper(SeasonList[[y]],Roter_par)
    
    Roter_PhenoFlex_predictions$Error<-
      Roter_PhenoFlex_predictions$predicted-Roter_PhenoFlex_predictions$pheno
    
    ggplot(Roter_PhenoFlex_predictions,aes(x=pheno,y=predicted)) +
      geom_point() +
      geom_abline(intercept=0, slope=1) +
      theme_bw(base_size = 15) +
      xlab("Observed bloom date (Day of the year)") +
      ylab("Predicted bloom date (Day of the year)") +
      ggtitle("Predicted vs. observed bloom dates")
    
    ggplot(Roter_PhenoFlex_predictions,aes(Error)) +
      geom_histogram() +
      ggtitle("Distribution of prediction errors")

  3. Compute the model performance metrics RMSEP, mean error and mean absolute error.

    library(chillR)
    library(ggplot2)
    
    
    Roter_par<-read_tab("data/PhenoFlex_parameters_Roter_Boskoop.csv")[,2]
    
    SeasonList <- genSeasonList(hourtemps$hourtemps, mrange = c(8, 6), years=c(1959:2018))
    
    Roter_PhenoFlex_predictions<-Roter_first[which(Roter_first$Year>1958),]
    
    for(y in 1:length(Roter_PhenoFlex_predictions$Year))
       Roter_PhenoFlex_predictions$predicted[y]<-PhenoFlex_GDHwrapper(SeasonList[[y]],Roter_par)
    
    Roter_PhenoFlex_predictions$Error<-
      Roter_PhenoFlex_predictions$predicted-Roter_PhenoFlex_predictions$pheno
    
    RMSEP(Roter_PhenoFlex_predictions$predicted,Roter_PhenoFlex_predictions$pheno)
    ##[1] 3.869168
    
    RPIQ(Roter_PhenoFlex_predictions$predicted,Roter_PhenoFlex_predictions$pheno)
    ##[1] 3.682963
    
    mean(Roter_PhenoFlex_predictions$Error)
    ##[1] -1.077083
    
    mean(abs(Roter_PhenoFlex_predictions$Error))
    ##[1] 2.859028

    RMSEP is 3.87, mean error is -1.08, and mean absolute error is 2.86.

28 The PhenoFlex model - a second look

Exercise:

Make chill and heat response plots for the ‘Alexander Lucas’ PhenoFlex model for the location you did the earlier analyses for.

Since the Phenoflex functions take too much time to load, I did not run the previous code example to get the phenology parameters for ‘Alexander Lucas’ dataset. I merely downloaded from here. I skipped the chill and heat response curves for ‘Alexander Lucas’ dataset because they are already presented in the module.

library(ggplot2)
library(tidyr)
require(reshape2)
library(chillR)

apply_const_temp <- function(temp, A0, A1, E0, E1, Tf, slope, portions=1200, deg_celsius=TRUE)
    {
  temp_vector <- rep(temp, times=portions)
  res <- chillR::DynModel_driver(temp=temp_vector,
                         A0=A0, A1=A1,
                         E0=E0, E1=E1,
                         Tf=Tf,
                         slope=slope,
                         deg_celsius=deg_celsius)
  return(res$y[length(res$y)])
}

gen_bell <- function(par, temp_values=seq(-5, 20, 0.1)) {
  E0 <- par[5]
  E1 <- par[6]
  A0 <- par[7]
  A1 <- par[8]
  Tf <- par[9]
  slope <- par[12]

  y <- c()
  for(i in seq_along(temp_values)) {
    y[i] <- apply_const_temp(temp=temp_values[i],
                             A0=A0, A1=A1, E0=E0, E1=E1, Tf=Tf, slope=slope)
  }
  return(invisible(y))
}

GDH_response<-function(T, par)
{Tb<-par[11]
 Tu<-par[4]
 Tc<-par[10]
 GDH_weight <- rep(0, length(T))
 GDH_weight[which(T >= Tb & T <= Tu)] <-
   1/2 * (1 + cos(pi + pi * (T[which(T >= Tb & T <= Tu)] - Tb)/(Tu - Tb)))
 GDH_weight[which(T > Tu & T <= Tc)] <-
   (1 + cos(pi/2 + pi/2 * (T[which(T >  Tu & T <= Tc)] -Tu)/(Tc - Tu)))
  return(GDH_weight)
}



Alex_par<-read_tab("data/PhenoFlex_parameters_Alexander_Lucas.csv")[,2]

temp_values=seq(-5, 30, 0.1)

temp_response<-data.frame(Temperature=temp_values,
                          Chill_response=gen_bell(Alex_par, temp_values),
                          Heat_response=GDH_response(temp_values,Alex_par))

melted_response<-melt(temp_response,id.vars="Temperature")


##for chill and heat efficiency
latitude<-48.5

month_range<-c(10,11,12,1,2,3)

Tmins=c(-20:20)
Tmaxs=c(-15:30)

mins<-NA
maxs<-NA
chill_eff<-NA
heat_eff<-NA
month<-NA

simulation_par<-Alex_par

library(lubridate)

for(mon in month_range)
    {
    weather<-make_all_day_table(data.frame(Year=c(2001,2001),
                                         Month=c(mon,mon),
                                         Day=c(1,days_in_month(mon)),
                                         Tmin=c(0,0),
                                         Tmax=c(0,0)))
    for(tmin in Tmins)
      for(tmax in Tmaxs)
        if(tmax>=tmin)
          {
          weather$Tmin<-tmin
          weather$Tmax<-tmax
          hourtemps<-stack_hourly_temps(weather,
                                        latitude=latitude)$hourtemps$Temp
          chill_eff<-c(chill_eff,
                       tail(PhenoFlex(temp=hourtemps,
                                 times=c(1: length(hourtemps)),
                                 A0=simulation_par[7],
                                 A1=simulation_par[8],
                                 E0=simulation_par[5],
                                 E1=simulation_par[6],
                                 Tf=simulation_par[9],
                                 slope=simulation_par[12],
                                 deg_celsius=TRUE,
                                 basic_output=FALSE)$y,1)/
                         days_in_month(mon))
                                       
          heat_eff<-c(heat_eff,
                      tail(cumsum(GDH_response(hourtemps,
                                          simulation_par)),1)/
                        days_in_month(mon))
          mins<-c(mins,tmin)
          maxs<-c(maxs,tmax)
          month<-c(month,mon)
        }
    }
results<-data.frame(Month=month,Tmin=mins,Tmax=maxs,Chill_eff=chill_eff,Heat_eff=heat_eff)
results<-results[!is.na(results$Month),]


write.csv(results,"data/model_sensitivity_PhenoFlex.csv")



Chill_sensitivity_temps<-function(chill_model_sensitivity_table,
                                  temperatures,
                                  temp_model,
                                  month_range=c(10,11,12,1,2,3),
                                  Tmins=c(-20:20),
                                  Tmaxs=c(-15:30),
                                  legend_label="Chill/day (CP)")
{
  library(ggplot2)
  library(colorRamps)

  cmst<-chill_model_sensitivity_table
  cmst<-cmst[which(cmst$Month %in% month_range),]
  cmst$Month_names<- factor(cmst$Month, levels=month_range,
                            labels=month.name[month_range])  
  
  DM_sensitivity<-ggplot(cmst,aes(x=Tmin,y=Tmax,fill=.data[[temp_model]])) +
    geom_tile() +
    scale_fill_gradientn(colours=alpha(matlab.like(15), alpha = .5),
                         name=legend_label) +
    xlim(Tmins[1],Tmins[length(Tmins)]) +
    ylim(Tmaxs[1],Tmaxs[length(Tmaxs)])
  
  temperatures<-
    temperatures[which(temperatures$Month %in% month_range),]
  temperatures[which(temperatures$Tmax<temperatures$Tmin),
               c("Tmax","Tmin")]<-NA
  temperatures$Month_names <- factor(temperatures$Month,
                                     levels=month_range, labels=month.name[month_range])  
  
  DM_sensitivity +
    geom_point(data=temperatures,
               aes(x=Tmin,y=Tmax,fill=NULL,color="Temperature"),
               size=0.2) +
    facet_wrap(vars(Month_names)) +
    scale_color_manual(values = "black",
                       labels = "Daily temperature \nextremes (°C)",
                       name="Observed at site" ) +
    guides(fill = guide_colorbar(order = 1),
           color = guide_legend(order = 2)) +
    ylab("Tmax (°C)") +
    xlab("Tmin (°C)") + 
    theme_bw(base_size=15)

}


Model_sensitivities_PhenoFlex<-
  read.csv("data/model_sensitivity_PhenoFlex.csv")
MApfelstetten_weather<-read_tab("data/MApfelstetten_weather.csv")


Chill_sensitivity_temps(Model_sensitivities_PhenoFlex,
                        MApfelstetten_weather,
                        temp_model="Chill_eff",
                        month_range=c(10,11,12,1,2,3),
                        Tmins=c(-20:20),
                        Tmaxs=c(-15:30),
                        legend_label="Chill per day \n(arbitrary)") +
  ggtitle("PhenoFlex chill efficiency ('Alexander Lucas')")




Chill_sensitivity_temps(Model_sensitivities_PhenoFlex,
                        MApfelstetten_weather,
                        temp_model="Heat_eff",
                        month_range=c(10,11,12,1,2,3),
                        Tmins=c(-20:20),
                        Tmaxs=c(-15:30),
                        legend_label="Heat per day \n(arbitrary)") +
  ggtitle("PhenoFlex heat efficiency ('Alexander Lucas')")

29 Can we improve the performance of PhenoFlex?

Exercise:

  1. What was the objective of this work?

    To test and recalibrate parameters of Phenoflex under different conditions.

  2. What was the main conclusion?

    Phenoflex, like every other model has limitations; it may be that under extreme conditions (very cold and hot temperatures), the model fails to include the necessary parameters taking play in breaking dormancy. More tests are needed.

  3. What experiments could we conduct to test the hypothesis that emerged at the end of the conclusion?

    Phenoflex should be tested more with observed data from very cold, very hot climates, and in between (progression), to be able to clearly draw the validity domains.

30 Frost risk analysis

Exercise:

  1. Download the phenology dataset for the apple cultivar Roter Boskoop from Klein-Altendorf.

    V1

    Pheno_year,First_bloom,Full_bloom,Last_bloom

    1958,19580509,19580511,19580514

    1959,19590414,19590421,19590504

    1960,19600424,19600505,19600510

    1961,19610409,19610419,19610501

    1962,19620514,19620522,19620530

  2. Illustrate the development of the bloom period over the duration of the weather record. Use multiple ways to show this - feel free to be creative.

    library(chillR)
    library(leaflet)
    library(dplyr)
    library(reshape2)
    library(kableExtra)
    library(ggplot2)
    library(Kendall)
    
    
    CKA_Roter<-read.csv("data/Roter_Boskoop_bloom_1958_2019.csv")
    CKA_weather<-read.csv("data/TMaxTMin1958-2019_patched.csv")
    
    Roter<-melt(CKA_Roter,
                          id.vars = "Pheno_year",
                          value.name="YEARMODA")
    Roter$Year<-as.numeric(substr(Roter$YEARMODA,1,4))
    Roter$Month<-as.numeric(substr(Roter$YEARMODA,5,6))
    Roter$Day<-as.numeric(substr(Roter$YEARMODA,7,8))
    Roter<-make_JDay(Roter)
    
    ggplot(data=Roter,aes(Pheno_year,JDay,col=variable)) +
      geom_line() +
      theme_bw(base_size=15) +
      scale_color_discrete(
        name="Phenological event",
        labels=c("First bloom", "Full bloom", "Last bloom")) +
      xlab("Phenological year") +
      ylab("Julian date (day of the year)")
    
    
    ggplot(data=Roter,aes(Pheno_year,JDay,col=variable)) +
      geom_line() +
      theme_bw(base_size=15) +
      scale_color_discrete(name="Phenological event", labels=c("First bloom", "Full bloom", "Last bloom")) +
      xlab("Phenological year") +
      ylab("Julian date (day of the year)") +
      geom_smooth(method = "lm")
    
    ggplot(data=Roter,aes(Pheno_year,JDay,col=variable)) +
      geom_smooth() +
      theme_bw(base_size=15) +
      scale_color_discrete(
        name="Phenological event",
        labels=c("First bloom", "Full bloom", "Last bloom")) +
      xlab("Phenological year") +
      ylab("Julian date (day of the year)") 

  3. Evaluate the occurrence of frost events at Klein-Altendorf since 1958. Illustrate this in a plot.

    library(chillR)
    library(leaflet)
    library(dplyr)
    library(reshape2)
    library(kableExtra)
    library(ggplot2)
    library(Kendall)
    
    
    CKA_weather<-read.csv("data/TMaxTMin1958-2019_patched.csv")
    
    frost_df=data.frame(
      lower=c(-1000,0),
      upper=c(0,1000),
      weight=c(1,0))
    
    frost_model<-function(x) step_model(x,frost_df)
    
    hourly<-stack_hourly_temps(CKA_weather,latitude=50.625)
    frost<-tempResponse(hourly,models=c(frost=frost_model))
    
    frost_model_no_summ<-function(x) step_model(x, frost_df, summ=FALSE)
    
    hourly$hourtemps[,"frost"]<-frost_model_no_summ(hourly$hourtemps$Temp)
    
    Daily_frost_hours<-aggregate(hourly$hourtemps$frost,
                                 by=list(hourly$hourtemps$YEARMODA),
                                 FUN=sum)
    
    Daily_frost<-make_JDay(CKA_weather)
    
    Daily_frost[,"Frost_hours"]<-Daily_frost_hours$x
    
    
    Daily_frost$Frost_hours[which(Daily_frost$Frost_hours==0)]<-NA
    
    ggplot(data=Daily_frost,aes(Year,JDay,size=Frost_hours)) +
      geom_point(col="light blue",alpha = 0.8) + 
      scale_size(range = c(0, 3),
                 breaks = c(1, 5, 10, 15, 20),
                 labels = c("1", "5", "10", "15", "20"),
                 name="Frost hours") +
      theme_bw(base_size=15)

  4. Produce an illustration of the relationship between spring frost events and the bloom period of ‘Roter Boskoop.’

    library(chillR)
    library(leaflet)
    library(dplyr)
    library(reshape2)
    library(kableExtra)
    library(ggplot2)
    library(Kendall)
    
    
    CKA_Roter<-read.csv("data/Roter_Boskoop_bloom_1958_2019.csv")
    CKA_weather<-read.csv("data/TMaxTMin1958-2019_patched.csv")
    
    #bloom
    Roter<-melt(CKA_Roter,
                          id.vars = "Pheno_year",
                          value.name="YEARMODA")
    Roter$Year<-as.numeric(substr(Roter$YEARMODA,1,4))
    Roter$Month<-as.numeric(substr(Roter$YEARMODA,5,6))
    Roter$Day<-as.numeric(substr(Roter$YEARMODA,7,8))
    Roter<-make_JDay(Roter)
    
    ##frost
    frost_df=data.frame(
      lower=c(-1000,0),
      upper=c(0,1000),
      weight=c(1,0))
    
    frost_model<-function(x) step_model(x,frost_df)
    
    hourly<-stack_hourly_temps(CKA_weather,latitude=50.625)
    frost<-tempResponse(hourly,models=c(frost=frost_model))
    
    frost_model_no_summ<-function(x) step_model(x, frost_df, summ=FALSE)
    
    hourly$hourtemps[,"frost"]<-frost_model_no_summ(hourly$hourtemps$Temp)
    
    Daily_frost_hours<-aggregate(hourly$hourtemps$frost,
                                 by=list(hourly$hourtemps$YEARMODA),
                                 FUN=sum)
    
    Daily_frost<-make_JDay(CKA_weather)
    
    Daily_frost[,"Frost_hours"]<-Daily_frost_hours$x
    
    
    Daily_frost$Frost_hours[which(Daily_frost$Frost_hours==0)]<-NA
    
    ggplot(data=Roter,
           aes(Pheno_year,JDay,col=variable)) +
      geom_line() +
      theme_bw(base_size=15) +
      scale_color_discrete(
        name="Phenological event",
        labels=c("First bloom", "Full bloom", "Last bloom")) +
      xlab("Phenological year") +
      ylab("Julian date (day of the year)") +
      geom_point(data=Daily_frost,
                 aes(Year,JDay,size=Frost_hours),
                 col="light blue",alpha = 0.8) + 
      scale_size(range = c(0, 3),
                 breaks = c(1, 5, 10, 15, 20),
                 labels = c("1", "5", "10", "15", "20"),
                 name="Frost hours") +
      theme_bw(base_size=15)
    
    
    
    ggplot(data=Roter,
           aes(Pheno_year,JDay,col=variable)) +
      geom_line() +
      theme_bw(base_size=15) +
      scale_color_discrete(
        name="Phenological event",
        labels=c("First bloom", "Full bloom", "Last bloom")) +
      xlab("Phenological year") +
      ylab("Julian date (day of the year)") +
      geom_point(data=Daily_frost,
                 aes(Year,JDay,size=Frost_hours),
                 col="light blue",
                 alpha = 0.8) + 
      scale_size(range = c(0, 3),
                 breaks = c(1, 5, 10, 15, 20),
                 labels = c("1", "5", "10", "15", "20"),
                 name="Frost hours") +
      theme_bw(base_size=15) +
      ylim(c(77,140))
    
    
    Ribbon_Roter<-dcast(
      Roter,Pheno_year ~ variable, value.var = "JDay")
    
    ggplot(data=Ribbon_Roter,aes(Pheno_year)) +
      geom_ribbon(aes(ymin = First_bloom, ymax = Last_bloom),
                  fill = "light gray") +
      geom_line(aes(y = Full_bloom)) +
      theme_bw(base_size=15) +
      xlab("Phenological year") +
      ylab("Julian date (day of the year)") +
      geom_point(data=Daily_frost,
                 aes(Year,JDay,size=Frost_hours),
                 col="light blue",
                 alpha = 0.8) + 
      scale_size(range = c(0, 3),
                 breaks = c(1, 5, 10, 15, 20),
                 labels = c("1", "5", "10", "15", "20"),
                 name="Frost hours") +
      theme_bw(base_size=15) +
      ylim(c(77,140))
    
    # identify frost events that overlap with bloom
    lookup_dates<-Ribbon_Roter
    row.names(lookup_dates)<-lookup_dates$Pheno_year
    
    Daily_frost[,"First_bloom"]<-
      lookup_dates[as.character(Daily_frost$Year),"First_bloom"]
    Daily_frost[,"Last_bloom"]<-
      lookup_dates[as.character(Daily_frost$Year),"Last_bloom"]
    
    Daily_frost[
      which(!is.na(Daily_frost$Frost_hours)),"Bloom_frost"]<-
      "Before bloom"
    Daily_frost[
      which(Daily_frost$JDay>=Daily_frost$First_bloom),"Bloom_frost"]<-
      "During bloom"
    Daily_frost[
      which(Daily_frost$JDay>Daily_frost$Last_bloom),"Bloom_frost"]<-
      "After bloom"
    Daily_frost[
      which(Daily_frost$JDay>180),"Bloom_frost"]<-
      "Before bloom"
    
    ggplot(data=Ribbon_Roter,aes(Pheno_year)) +
      geom_ribbon(aes(ymin = First_bloom, ymax = Last_bloom),
                  fill = "light gray") +
      geom_line(aes(y = Full_bloom)) +
      theme_bw(base_size=15) +
      xlab("Phenological year") +
      ylab("Julian date (day of the year)") +
      geom_point(data=Daily_frost,
                 aes(Year,JDay,size=Frost_hours,col=Bloom_frost),
                 alpha = 0.8) + 
      scale_size(range = c(0, 5),
                 breaks = c(1, 5, 10, 15, 20),
                 labels = c("1", "5", "10", "15", "20"),
                 name="Frost hours") +
      scale_color_manual(
        breaks=c("Before bloom", "During bloom", "After bloom"),
        values=c("light green","red","light blue"),
        name="Frost timing") +
      theme_bw(base_size=15) +
      ylim(c(77,140))

  5. Evaluate how the risk of spring frost for this cultivar has changed over time. Has there been a significant trend?

    library(chillR)
    library(leaflet)
    library(dplyr)
    library(reshape2)
    library(kableExtra)
    library(ggplot2)
    library(Kendall)
    
    CKA_Roter<-read.csv("data/Roter_Boskoop_bloom_1958_2019.csv")
    CKA_weather<-read.csv("data/TMaxTMin1958-2019_patched.csv")
    
    #bloom
    Roter<-melt(CKA_Roter,
                          id.vars = "Pheno_year",
                          value.name="YEARMODA")
    Roter$Year<-as.numeric(substr(Roter$YEARMODA,1,4))
    Roter$Month<-as.numeric(substr(Roter$YEARMODA,5,6))
    Roter$Day<-as.numeric(substr(Roter$YEARMODA,7,8))
    Roter<-make_JDay(Roter)
    
    ##frost
    frost_df=data.frame(
      lower=c(-1000,0),
      upper=c(0,1000),
      weight=c(1,0))
    
    frost_model<-function(x) step_model(x,frost_df)
    
    hourly<-stack_hourly_temps(CKA_weather,latitude=50.625)
    frost<-tempResponse(hourly,models=c(frost=frost_model))
    
    frost_model_no_summ<-function(x) step_model(x, frost_df, summ=FALSE)
    
    hourly$hourtemps[,"frost"]<-frost_model_no_summ(hourly$hourtemps$Temp)
    
    Daily_frost_hours<-aggregate(hourly$hourtemps$frost,
                                 by=list(hourly$hourtemps$YEARMODA),
                                 FUN=sum)
    
    Daily_frost<-make_JDay(CKA_weather)
    
    Daily_frost[,"Frost_hours"]<-Daily_frost_hours$x
    
    
    Daily_frost$Frost_hours[which(Daily_frost$Frost_hours==0)]<-NA
    
    # identify frost events that overlap with bloom
    lookup_dates<-Ribbon_Roter
    row.names(lookup_dates)<-lookup_dates$Pheno_year
    
    Daily_frost[,"First_bloom"]<-
      lookup_dates[as.character(Daily_frost$Year),"First_bloom"]
    Daily_frost[,"Last_bloom"]<-
      lookup_dates[as.character(Daily_frost$Year),"Last_bloom"]
    
    Daily_frost[
      which(!is.na(Daily_frost$Frost_hours)),"Bloom_frost"]<-
      "Before bloom"
    Daily_frost[
      which(Daily_frost$JDay>=Daily_frost$First_bloom),"Bloom_frost"]<-
      "During bloom"
    Daily_frost[
      which(Daily_frost$JDay>Daily_frost$Last_bloom),"Bloom_frost"]<-
      "After bloom"
    Daily_frost[
      which(Daily_frost$JDay>180),"Bloom_frost"]<-
      "Before bloom"
    
    ##checking for bloom trends
    Bloom_frost_trend<-aggregate(
      Daily_frost$Frost_hours,
      by=list(Daily_frost$Year,Daily_frost$Bloom_frost),
      FUN=function(x) sum(x,na.rm=TRUE))
    colnames(Bloom_frost_trend)<-c("Year","Frost_timing","Frost_hours")
    
    DuringBloom<-
      Bloom_frost_trend[which(Bloom_frost_trend$Frost_timing=="During bloom"),]
    
    ggplot(data=DuringBloom,aes(Year,Frost_hours)) +
      geom_col() +
      ylab("Frost hours")
    
    
    Kendall(x=DuringBloom$Year,y=DuringBloom$Frost_hours)
    ##tau = 0.0333, 2-sided pvalue =0.74038
    
    lm(DuringBloom$Frost_hours~DuringBloom$Year)
    ##Call:
    #lm(formula = DuringBloom$Frost_hours ~ DuringBloom$Year)
    
    #Coefficients:
       #  (Intercept)  DuringBloom$Year  
     #      -91.48652           0.04786 

    Trend above is hard to interpret, I will not attempt to. And after checking via Kendall test, trend is insignificant (p-value >0.05) anyway.

Last words (for now)

I am glad I enlisted in this course. I learned from bottom and now I can say I am still at the bottom but a bit with clue what is going on (in terms of R, climate scenarios, and tree phenology). The errors I encountered forced me to dissect function per function, objects, lists, files, and the code-- it was dreadful but necessary.
Upon using ‘chillR’, I observed some things and I list them down below, hoping the developers will take a look at them for future improvements:

#'Warning: scenario_type doesn't look right - is this really a relative scenario?Warning: updating scenario_type to 'absolute'Error: temperature_scenario is NULL

With every word, I felt the passion that the module creators had while crafting the modules. It was very well-made to the point that someone clueless like I exited the course with huge progress in terms of R knowledge. I want more exposure to these concepts so I will be seeing the team again next semester for the Decisions Analysis and Forecasting subject. I will strive to do better for sure. Bis bald!

References

Jaffrés, Jasmine. (2019). GHCN-Daily: a treasure trove of climate data awaiting discovery. Computers & Geosciences. 122. 35-44. 10.1016/j.cageo.2018.07.003.

Glaser, Rüdiger & Himmelsbach, Iso & Boesmeier, Annette. (2017). Climate of migration? How climate triggered migration from southwest Germany to North America during the 19th century. Climate of the Past. 13. 1573-1592. 10.5194/cp-13-1573-2017.

Kaufmann, H. and Blanke, M.M. (2019). Chilling requirements of Mediterranean fruit crops in a changing climate. Acta Hortic. 1242, 275-280. DOI: 10.17660/ActaHortic.2019.1242.38
https://doi.org/10.17660/ActaHortic.2019.1242.38

Luedeling, Eike, Evan H Girvetz, Mikhail A Semenov, and Patrick H Brown. 2011. "Climate Change Affects Winter Chill for Temperate Fruit and Nut Trees." PloS One 6 (5): e20155. https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0020155.

Milech, Chaiane & Scariotto, Silvia & Dini, Maximiliano & Herter, Flavio & Raseira, Maria. (2018). MODELS TO ESTIMATE CHILLING ACCUMULATION UNDER SUBTROPICAL CLIMATIC CONDITIONS IN BRAZIL. Revista Brasileira de Climatologia. 23. 10.5380/abclima.v23i0.53086. Tominaga, A., Ito, A., Sugiura, T., & Yamane, H. (2022). How Is Global Warming Affecting Fruit Tree Blooming? “Flowering (Dormancy) Disorder” in Japanese Pear (Pyrus pyrifolia) as a Case Study. Frontiers in Plant Science, 12. https://doi.org/10.3389/fpls.2021.787638

The Norns (2018) Norse Mythology for Smart People. Available at: https://norse-mythology.org/gods-and-creatures/others/the-norns/#:~:text=In%20Norse%20mythology%2C%20the%20Norns,any%20and%20all%20other%20beings. (Accessed: March 15, 2023).

Tominaga, A. et al. (2022) "How is global warming affecting Fruit Tree Blooming? 'Flowering (dormancy) disorder' in Japanese pear (Pyrus pyrifolia) as a case study," Frontiers in Plant Science, 12. Available at: https://doi.org/10.3389/fpls.2021.787638.

Schwäbisches Streuobstparadies (no date) Streuobstparadies. Available at: https://www.streuobstparadies.de/ (Accessed: March 15, 2023).

Xu, Yunjia & Dai, Junhu & Ge, Quansheng & Wang, Huanjiong & Tao, Zexing. (2020). Comparison of chilling and heat requirements for leaf unfolding in deciduous woody species in temperate and subtropical China. International Journal of Biometeorology. 10.1007/s00484-020-02007-7.